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ABSTRACT 


The electrical conductivity of the upper mantle is estimated 
from low latitude magnetic field variations (magnetic storms) caused 
by large fluctuations in the equatorial ring current. The data base 
is derived from satellite magnetic field measurements which offer 
better global coverage than land based observatories. 

The procedures of analysis consist of 1) separation of the 
disturbance field into internal and external parts relative to the 
surface of the earth, 11) estimation of a response function "Q(u))" 
which relates the internally generated magnetic field variations to 
the external variations due jto the ring current, and ill) interpretation 
of the estimated response function using theoretical response 
functions for known conductivity profiles. Special consideration is 
given to possible ocean effects. 

Magnetic field variations are derived from magnetic field 
magnitude data measured by satellites Ogo 2,4, and 6 which collected 
data fr jm October, 1965 to July, 1971. Flying nearly polar orbits, 
these satellites sampled the field every H minute over an altitude 
range of 400 km to 1500 km (Cain and Langel,1971; Langel, 1974). 

Best estimates of the response function Q(o)) were obtained 
using stacked, smoothed power spectra from eight storms. The frequency 
range for these smoothed estimates is from 0.2 cycles/day to 
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2.0 cycles/day. Also, to examine r possitiiy source of noise In the 
data, a model of the crustal anomaly field (model by Mayhew, 1980) 
was removed from the magnetic field data and a "corrected" Q was 
calculated. Although the corrected Q differed slightly from the 
original Q for individual storms, the stacked corrected Q did not 
differ from the uncorrected stacked results. It is concluded that within 
the reliability limits of the data, the anomaly field has no significant 
contribution. 

Using a finite difference algorithm, a theoretical Q(oj) 

is calculated for a given conductivity profile. Matching a theoretical 

Q to the measured Q implies a model where most of the upper mantle 

has a conductivity of order 10 mho/m. This value agrees with Banks' 

(1972) model but differs from Parker's (1970) value of order 10“^mho/m. 

Considering the high temperatures implied by Parker's model, geochemical 

and petrological evidence favors the 10 mho/m mantle. 

Since the frequency range (.2 to 2. cpd) extends beyond 

that used by Banks and Parker, considerations were made for possible 

ocean effects. Representing the top 3 km of the earth with a shell of 

conductivity equal to a weighted sum of oceanic and continental 

conductivities, the ocean effect did not appear in our frequency 

_2 

range until the continental conductivity was raised to 10 mho/m — this 

_2 

sets an upper limit of 10 mho/m on the top 3 km of crust. Also suggested 

—3 

by the data is a lower crust-upper mantle conductivity of 10 raho/m , 
extending downward to 30 km (upper limit) . 


A temperature profile is proposed using conductivity*- 


temperature data for single crystal olivine ('ViFo^q) (Duba» J976). 

The resulting temperature profile is reasonable for depths below 150- 
200 km, but is too high for shallower depths. Apparently, conductivity 
is not controlled solely by olivine at shallow depths. 

As a new application of satellite data, the results are 
most promising. Reliability of the estimates of Q(u) will be greatly 
Improved if more data sets are stacked. 
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CHAPTER Is INTRODUCTION 


This thctls concerns the earth's electrical conductivity as 
determined from low latitude magnetic field variations (magnetic 
storms) caused by large fluctuations in the equitorial ring current. 
Using satellite observations of these magnetic field variations, 
information about the electrical conductivity structure and the 
associated temperature profile can be obtained. 

The procedures of analysis consist of 1) separation of the 
disturbance field into internal and external parts relative to the 
surface of the earth, ii) estimation of a response function "Q(w)" 
which relates the internally generated magnetic field variations to 
the external variations due to the ring current, and ill) interpretation 
of the estimated response function using theoretical response functions 
for known conductivity profiles. Special consideration is given to 
possible ocean effects. 

I.l Previous investigations. 

In 1889, Schuster used observations of the diurnal (Sq) 
magnetic variations to infer the general conductivity structure of the 
earth. Using Gauss's spherical harmonic representation of magnetic 
scales potential 

U - a E [1® + e® (|)^] P®(e,(f!) (1.1) 

il,m 

where a is the radius of the earth, r is the radial distance of the 
observer, i® is the internal source coefficient and e® is the external 
source coefficient for spherical harmonic P®(0,4>) of degree I and 
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ord.r m, Sohuit.r .•clm.t.d .nd and coiclud.d that tha aarth 

could be re^: .eeentcd by an equivalent conducting iphere with a amaller 

radius than the Earth's. Chapman (19 L9) improved upon Schuster's 

analysis and concluded that the equivalent sphere has a conductivity 
•2 

of 3. X 10 mho/m and a radius of 250 km less than that of the earth. 

Chapman and Price (1930) also studied diurnal variations and 

obtained estimates of electrical conductivity compatible with Chapman's 

(1919) earlier estimates. However, in their study of a non-periodic 

phenomenon, namely magnetic storms, they found a significantly higher 

-1 

conductivity estimate of ah:?ut 4. x 10 mho/m. This higher conductivity 
estimate was taken as an indication of a non-uniform conductivity 
structure. 

Lahlri and Price (1939) further developed the theory for non- 
uniform conductors and obtained results which supported Chapman and 
Price's view that an Increase in conductivity existed at a depth of 
250 km. They also concluded that an additional large increase in 
conductivity (one to two orders of magnitude) exists at a depth of 
700 km. Their results indicated an enhanced conductivity near the 
surface of the earth, which they attributed to the oceans. 

To Improve the resolution of the electrical conductivity 
structure, analysis of more magnetic variations of differing fre- 
quencies was needed. Intuitively, this becomes apparent when one 

k 

considers the electromagnetic skin depth, 6 ■ (2/woy) , which is 
inversely proportional to frequency oi, conductivity a and magnetic 
permeability y. To obtain resolution near the surface of the earth. 


orui nttdi magnetic variations with small skin depths. These have high 
frequency components which are attenuated before reaching far into 
the Interior, To probe deeper into the earth • the magnetic variations 
must have larger skin depths and therefore lower frequency components 
so they can reach the interior before they are greatly attenuated. 
Theoretically, a necessary condition for obtaining a unique model of 
the earth's conductivity structure is that magnetic variations with 
an Infinite range of frequencies must be sampled (Bailey, 1970). 

However , an infinite frequency range of magnetic variations is 

not available for Induction studies. At the low frequency end, below 
-3 

10 cpd (cycles per day) the geomagnetic spectrum is dominated by 
"secular" variations originating internally at the core-mantle boundary, 
At the high frequency end, above 0.25 cpd, lack of global coverage 
along with possible lateral variations in the upper 400 km of the 
earth hinders application of induction theory (Anderson, et al,, 1979). 

Not until the application of cross-spectral analysis tech- 
niques (Banks, 1969, 1972) could the Intervening continuum as well 
as some major peaks associated with periodic phenomena (l.e. "Sq") 
be included in induction studies. Using Gauss's potential field 
representation (eqn. 1.1), Banks separated the observed magnetic 
variations into internal and external origin. He then defined a 
frequency response function Q™(w) by 
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where l”(w) and E™(o)) are the frequency spectra for the Internal and 
external coefficients of degree i and order m In eqn. 1.1. Using 
standard cross^spectral techniques (Blackman and Tukey, 1958), Banks 
estimated the response function (mainly Q^(w) harmonic) and conducted 
a Monte Carlo search fo* the most compatible conducting profile 
(fig. 1.1). 

Banks' estimate of Q^(w) and his resulting conductivity profile 
are based on his analysis of magnetic variations of frequencies 
0.003 cpd to 0.25 cpd. For frequencies above 0.25, estimates of Q^(o)) 
had nonphysical, negative phases. He attributed thia lack of relia- 
bility for high frequency estimates to lack of global coverage and/or 
lateral variations In the upper mantle and crust. 

Using Banks' data, Parker (1970) suggested another conductivity 
profile. Pa. xer obtained this profile (fig. 1.2) by applying Inverse 
theory after the fashion of Backus and Gilbert. As compared to Banks' 
original 1969 model, Parker's model confirms the sharp rise in con- 
ductivity at 400 km but differs In its near surface value: near 0.1 

—2 

mho/m as opposed to Banks' figure of 10 mho/m. 

In the models thus far presented, existence of the oceans is 
not considered. Since the high frequency estimates of Q^(w), near 
0.25 cpd. Is from the zonal (local time independent) part of magnetic 
storm variations, 'Dst,' which is assumed to have only pj component 
(eqn. 1.1), the effect of oceans are considered unimportant for Dst 
variations (Lahirl and Price, 1939). For Sq variations (1 cpd and 
harmonics) which have higher harmonic components, the oceans probably 


Figure 1.1: Banks* 1972 conductivity profile. (J. Geomag, 


Geoelectr., 24, 337-351, 1972) 


Figure 1.2; Parker's 1970 conductivity profile. 


(Geophys . J . 


R. astr. Soc. (1970) 22, 121-138) 
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have a significant effect (Lahlrl and Price, 1939; Bullard and Parker, 
1970; Banks, 1972, Jady, 1974a, b). Banks (1972) later reworked his 
data to Include Sq variations, but he does not Include a layer for 
ocean conductivity. Ocean effects are usually Included in local studies 
only. More discussion of ocean efi.*,cts will be presented in chapter V. 

As noted earlier, the low frequency end of the geomagnetic 
spectrum is dominated by "secular variations." Since these variations 
are of purely Internal origin, induction methods used in the studies 
already described are not appropriate. However, secular variations 
still can yield information about the earth's deep interior. 

In a method developed by McDonald (1957) , the attenuation 
behavior of a sinusoidal magnetic variation originating at the core- 
mantle boundary is used as an indicator of lower mantle conductivity. 
Assuming spherical symmetry and a power law function for conducticlty, 
McDonald calculates ratios of attenuation functions of differing 
harmonic degree as would be observed at the earth's surface for differ- 
ent conductivity structures. Comparing these predicted ratios to 
those obtained from spherical harmonic representation of observed 

secular variations, McDonald inferred a conductivity profile with the 

2 

lower mantle conductivity of 10 mho/m. 

As McDonald's work estimates the conductivity of the lower 
mantle, local magnetotellunlc studies estimate crustal conductivities. 
Local measurements show the first few kilometers of the earth may be 
fairly conductive, typically 0.01 to 0.1 mho/m while the underlying 
rock may be as low as 10 ^ mho/m (Gough, 1973). Sea water varies in 
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conductivity according to temperature and salinity, but typically 
is 3.3 mho/m. 

The global coverage essential to conductivity studies of the 
upper 400 km of the earth can be provided by satellite observations. 
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CHAPTER II: SOURCES OF MAGNETIC VARIATIONS 

Part of the objective of this study is to estimate a response 
function (eqn, 1.2) relating the time-varying magnetic field external 
to a conductor to that generated internally. Chapman (1919) concluded 
that this "conduction" was within the earth and that the upper 250 km 
and the surrounding space were nonconductors. Later Investigations 
concluded that the upper 250 km did have significant conductivity 
and the boundary between conductor and nonconductor became the earth's 
surface. In this study, the boundary shall be the surface of the 
earth. This is a valid assumption if there are no electric currents 
(hence, no conductors) between the satellite and the surface of the 
earth. Ionospheric currents do exist between the satellite and the 
earth surface, but these will be demonstrated to be unimportant in 
this study. 

Another assumption in this study is that the equitorial ring 
current, represented by a P^ harmonic magnetic variation, is the sole 
source of external magnetic field variations. This implies that the 
observed magnetic variations of internal origin are due only to 
induction effects from the ring current. This assumption is valid if 
other current systems in the ionosphere and magnetosphere and their 
Induced fields are accounted for. 

To examine the validity of these assumptions, we shall consider 
the various electric current systems with regard to the position of 
the satellite. The trajectories of satellites Ogo-2, 4, and 6 place 
them in a polar orbit with altitudes varying 400 km to 1500 km. This ' 
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orbit situates the satellite abovre the denser regions of the ionos- 
phere and well Inside the magnetosphere. 

In the following section, electric current systems will be 
presented as "external" or "internal" to the satellite's orbit. 

II. 1 External current systems . 

The equatorial ring current and Its associated magnetic 
variations provides the source field for this study. Represented 
approximately by a toroidal current flow, Its Influence on magnetic 
field variations during magnetic storms is fairly uniform over the 
mid- and low latitudes of the earth. However, as the earth's magnetic 
field Is distorted towards the sunward side, the ring current Is also 
asymmetric In shape with respect to the sun. This causes a slight 
local time dependence for these variations. As for location, Its 
Inner boundary Is generally located at a distance greater than 3 earth 
radii from the earth, though occasionally it will extend Into the upper 
Ionosphere (fig. 2.1). In any case. It Is located well outside of the 
orbits of the satellites. 

From a long tradition In characterizing land observatory data, 
the magnetic disturbance field "D" is expressed in terms of an axially 
symmetric part "Dst" and a longitudinal or local time dependent part 
"DS." Dst or "storm time variation" is the average value of the change 
In the horizontal component of the disturbance field measured along 
the equator. This Is the principal Influence of the equitorlal ring 
current for low and mid-latitudes. DS or "disturbance local-time 
inequality" is the variation in D which depends on local time (Rishbeth 


Figure 2,1; (a) The earth's magnetosphere as viewed from the 
equatorial plane. 

(b) Location of the Ionosphere In the atmosphere. 

(Rishbeth and Garriott, Introduction to Ionospheric 
Physics , Academic Press, New York, 1969) 


Temperature T Electron 11 

Density N 1 0 rn 
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and Garrlott, 1969). 

D - Dst(t) + DS 

where (2.1) 

DS •> 2 ? sin (nA + s' ) 
n n n 

t is univetsal time, A is longitude (relative to the sun), and 
and are the series coefficients. For studies using land obser- 
vatory data, the DS component of the variation must be removed oefore 
analysis of the ring current variation can be made. 

However, with satellite observations, the trajectories are 
such that measurements over the mid- and low latitudes are at a 
constant local time. Hence, the disturbance field is measured at a 
constant local time over successive orbits (see Ch. Ill) and local 
time variations need not be removed. Unlike land observatories, the 
satellite is constantly changing in latitude "6", so different sources 
of magnetic variations are measured. This is particularly true in the 
polar regions of the earth. 

Current systems in the auroral zones are exceptionally active 
during magnetic disturbances. During such times, charged particles 
from the upper ionosphere and magnetosphere are absorbed in the polar 
regions of the earth. These charged particles contribute to the 
"polar electrojet" and other related current systems causing increased 
magnetic variations in the high latitudes. Current systems associated 
with auroral phenomena are more complicated in geometry and in tem- 
poral behavior than the equatorial ring current. They also involve 
regions of space both above and below the orbits of the satellites 
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and hanct art not aultabla for induction atudias, at laaat not in tha 
mannar praaantad hara. Tharaforai tha baat way to daal with polar 
currant syatama ia to axcluda magnatic data takan in tha highar 3ati- 
tudas. 

II. 2 Intarnal currents . 

The electric currants induced within the earth by tha magnatic 
variations due to the ring current are the currants of interest in 
this study. However, other currents exist in the earth and in the 
intervening ionosphere between the earth and the satellite. 

The daytime ionospheric current systems originate in the "E- 
layer" of the ionosphere at altitude ranges of 90-130 km. During the 
course of a day, the tidal effects of the sun (and moon) produce 
electric currents in the ionosphere. This "Sq" system causes the 
diurnal variations observed in dally magnetic records with its peak 
disturbance at local noon. Sq variations also Induce electric 
currents within the earth and oceans that contribute roughly one- 
third to the total observed magnetic variation attributed to Sq 
(Rishbeth and Garrlott, 1969). 

The geometry of the Sq current system is very complicated 
and, when given a spherical harmonic representation, requires several 
degrees of spherical harmonica. The currents induced in the earth, 
particularly in the oceans, also are complicated in shape. In their 
study of ocean effects, Bullard and Parker (1970) have analysed the 
effect of Sq variations induced within the oceans. Because of the 
ocean shapes and the insulating effects of the continents, the main 
influence of Sq is in harmonics higher than P^. 
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Sq Is assoclatsd not only with ths tidal affacts of tha sun, 
but also with tha danslty of fraa alactrons In tha lonosphara; tha 
dagraa of Ionisation and, tharafora, tha magnatlc variations dua to 
Sq ara local'*tlma dapandant (contrlbuta to DS In aquation 2.1), and 
ara at a mlnlmuin naar local midnight. 

Sometimes considered part of the Sq system Is the "equatorial 
electrojet." This current flows eastward along the dipole equator 
during the day, causing large magnetic variations at equatorial land 
observatories. A return current flows at night, but the affect Is 
much less Intense than during tie day. 

These "quiet day variations," Sq, L (lunar variations), and 
the equatorial electrojet, occur every day and are periodic in nature. 
They are easily detected by land observatories. Their effects on the 
data can bs minimized by applying a one-cycle/day filter to the data 
or by simply considering data at one local time, particularly at night. 

During magnetic disturbances, additional currents appear or 
are enhanced In the Ionosphere; In the polar regions, current systems 
associated with Increased auroral activity are stimulated. Higher 
electron densities appear in the E-layer, principally in the high 
latitudes. Other currents are Induced in lower latitudes, but these 
are strongly local-time dependent because of solar effects on ioniza- 
tion. It then follows that DS variations from Ionospheric currents 
are minimized during the night for the low and mid-latitudes; in the 
polar regions where particle absorption from the upper ionosphere and 
the magnetosphere takes place, magnetic variations can be active in 
nighttime as well as the daytime. 
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Th. pj m.r..tlc v.rl.tion. du. to th. ring cutr.nt Indue, 
currsnti within th* ionoiph«r« . How«v*r , th* int*gr«t*d conductlviti** 
for th* iono*ph*r*g meaturi^d ov*r hundr«d* of kiXom«t*r*t vary fros 10 
to 1,000 mhos (dapandant on diraction of intagration) and tharefora 
imply a vary low conductanca with ragard to this study. Induction 
offsets in th* ionosphara will ba assantially nois* in tha data. 

As for induction in th* ocaan from ring currant fluctuations, 

0 12 
th* currants should ba much lass than th* and currant* asso- 

ciatad with Sq. This raducad offset is bacause of tha insulating 

influanca of tha cont inants on a currant that would have to flow 

around the equitorial regions of the earth. The ocean effect is 

further discussed in Chapter V. 
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CHAPTER Ills DATA 

IIl.l Data collection . 

Tha magnltuda of tha magnatlc fiald was maasurad by satalllta^ 
Ogo-2, 4, and 6 which collectad data from Octohar, 1965 to July, 1971. 
Flying narrly polar orbits, thana tatallibas samplad tha fiald avary 
1/2 minuta ovar an altltuda ran{;a of 400 km to 1500 km (Cain and 
Langel, 1971; Langal, 1974). Theaa magnetic fiald maasuramants wars 
made by a rubidium vapor magnetometer with an accuracy of 6 gammas 
(lO"^ rrausa « lO”^ tesla.) . 

With each orbit taking approximately 90 minutes, the flight 
path from pole to pole essentially followed a longitudinal line. 

Since it took about 20 minutes to fly between -50* to +50® latitude, 
measurements at mid- and low latitudes were made at a constant local 
time (relative position of the sun). During the course of one orbit, 
the satellite crossed the equator twice and hence can be characterized 
by two local times corresponding to the two equitorlal crossings. 

Over a period of a few days, these local times for equatorial crossings 
remain almost constant. 

Data types for these sate! ,tes included geographic longitude 
and latitude, geomagnetic latitude, local time, universal tim';^, and 
altitude for each data point. This arrangement allows for easy 
elimination of data collected over the polar regions. It also allows 
for selection of data at particular local times. 
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111. 2 Data selection . 

In studying magnetic field variations associated with 
fluctuations In the equitorlal ring current, one needs data drawn 
from periods of large magnetic disturbances. Such data are here 
defined as having negative deviations In magnetic field magnitude 
greater than 50 gammas as compared to normal dally variations. 

Fu»'‘;her selection of data was based on local time. Since 
the satellite orbit Is characterized by two local times which remain 
constant over a few days, one has the option of selecting only that 
data which corresponds to one local time. Since the field fluctua- 
tions due to the ring current have a local time dependence, only 
data from one local time can be used. Furthermore, the nighttime 
data has been selected in order to minimize the effects of Ionization 
currents. 

The data and local times for data used in this study are 
presented in Table 3.1. 

111. 3 Field separation Into internal and external coefficients . 

For each data set listed in Tables 3.1 , two time 

series e(t ) and l(t ) are extracted from the satellite data. These 
n n 

time series correspond to the external coefficient "ej" and Internal 
coefficient "1^" of a harmonic representation of the scalar potential 
for disturbance field D 

6 « -vu 

U - a(ej (|) + ij (|)S pj (COS0) 


(3.1) 
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Date 

March 5» 1970 
July 26. 1969 
*April 21, 1970 
June 11. 1968 
December 31. 1967 
October 30. 1968 
February 8, 1968 
September 29. 1969 
tMarch 22. 1966 


TABLE 3.1: DATA 

Modified Julian Date 
40650-40657 
40428-40430 
40697-40698 
40018-40022 
39855-39860 
40159-40165 
39894-39900 
40493-40498 
39206-39209 


Local Time (Hrs.) 
19.50 
0.50 
1.75 

23.00 
15.75 

21.00 
23.00 
16.25 
15.75 


* included only in stacked data with anomaly correction, 
t included only in stacked date without anomaly correction. 
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where a Is the radius of the earth, r Is the radial distance of the 
observer, 9 is the colatitude and P® (cos 0 ) is the spherical harmonic 
of degree 1 and order 0. During the course of the magnetic distur- 
bance, the coefficients e® and i® change with universal time t, so 
the time series for e^ and 1^ can be defined. 

To identify these time series, the disturbance field D must 
be extracted from the observed field magnitude data |1°|. This is 
done by defining an observed change in total field, "AB°" by 

AB° - |B°1 - |b®| (3.2) 

where |b®{ corresponds to the quiet day field represented by a field 
model of degree 12 (Langel, 1974). A theoretical AB*^, which includes 
disturbance field 5 , 

AB*= - |b“ + d| - | 5 ®| (3.3) 

is then fitted to the observed AB°. Since D varies with colatitude 9 , 

this fit is performed for data points between -50® to 50* for each 

orbit. This results in one estimate of e^ and ij per orbit. When 

data is considered for the whole magnetic storm, each e^ and 1^ is 

estimated at a different universal time t, defined to be the time of 

crossing the geomagnetic equator. 

The resulting time series e(t ) and i(t ) are equally spaced 

n n 

in time by approximately 90 minutes and have a total period equal to 
the duration of the storm. 

This e(t ) and i(t ) analysis is a slightly different approach 
n n 

than that applied to land observatory data. Since observations on 
land vary with local time, the DS part of the disturbance field is 
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subtracted from the observations — only Dst Is used for magnetic storm 
analysis. In this satellite data, local time Is constant, so the 
analysis Includes more than just Dst. 

An example of a field separation for the strong magnetic 
disturbance of March 5-9, 1970 is presented In fig. 3.1. 

A possible source of error In this field separation Is the 
crustal anomaly field. To explore this effect, a field model (courtesy 
of Mike Mayhew, 1980) is subtracted from the total field measurements 
before the field separation is performed. This anomaly field proved 
to have small effect on large magnetic disturbances and on stacked 
data (see Ch. IV), but has some effect on field separation for weaker 
magnetic disturbances. 


22 


Figure 3.1: Field separation Into external (e(t)) and Internal <l(t)) 
parts for March 5» 1970 magnetic storm. 
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CHAPTER IV I DATA REDUCTION AND ANALYSIS 

The magnetic field data measured by the satellites are 

expressed in terms of two time series e(t ) and i(t.)» n 1, . . . , 

n n 

N for each magnetic disturbance observed. The external field varia- 
tion represented by e(t^) is the "input" and the internal field 
variation represented by i(t^) is the "output" for the linear system 

i(t) - q(t) * e(t) (4.1) 

where "*" denotes convolution and "q(t)" is the linear Impulse function. 
This operation can also be represented in the frequency domain by 

I(o») - Q(u)) E(o)) (4.2) 

where I(o>) and E(u)) are the Fourier transforms of l(t) and e(t) and 
Q(o)) is the "frequency response function." To ascertain the conduc- 
tivity structure of the earth, we first must estimate Q(o)). 

From eqn. (4.2), the frequency response function can be defined 
by 

0 ( 0 .) «- 3 ) 

If the above definition were applied to the original time series e(t^) 

and i(t ), n - 1, . . . , N, the resulting estimate of Q(u ), n « 1, 
n n 

. . . , N/2 would have the highest resolution that the data could 
provide, but the least reliability. 

To improve reliability of an estimate of Q(w), one needs to 
smooth the data, either in the time domain or the frequency domain. 
However if one were to smooth the frequency spectra I(o)) and E(w) and 
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than estinate 0(u) with eqn. (4.3), one would have a very biased 
estimate of Q(a>) (Bendet and Plersol, 1971). Instead, a better 
estimate of Q(u)) Involves smoothed estimates of the cross-power 
spectrum power spectrum 


Q(w) 




(4.4) 


knowing that the cross-power spectrum can be expressed as 

* 

Cj^2(w) - E (w) I(u) 


(4.5) 


where denotes complex conjugate, and the power spectrum C^^(oo) is 



(w) 


E (o)) E(u) 


(4.6) 


one can easily see that eqn. (4.4) is an equivalent definition i:o 
eqn. (4.3). The difference between the two definitions is in the 
amount of bias when one smooths in the frequency domain. Eqn. (4.4) 
will yield a less biased estimate than eqn. (4.3). 

To obtain a reliable estimate of Q(o)) for each data set listed 
in Table 3.1 , the following procedure as recommended by Bendat 

and Plersol (1971) was followed: 

1) ptaproMss ttaie series e(t^) end n - 1 N. 

2) calculate the Fourier transforms of e(t ) and l(t ) 

n n 

3) calculate the raw estimates of the cross-power spectrum 

A A 

Ci2(w) and power spectrum Cj^j^(w) 

A A 

4) smooth the spectra frequencies 

i - 1 M 

5) estimate Q(oi^), i - 1,.., M, M < N/2. 
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The preproceaalng In step (1) involves removing a linear trend, 

defined as the line which was tangent to the endpoints of the time 

series. The resulting time series starts and ends at the zero base 

line with the maximum well above zero. 

The Fourier transforms of e(t ) and i(t ) are routinely calcu- 

n n 

lated using the fast Fourier transform or FFT algorithm. The resulting 
spectra have half the number of original data points due to folding 
around the zero frequency axis. 

A 

In step (4) the raw estimates C, .(od ) and C.,(u ), n ■ 1, . . . , 

n JLX n 

N/2 are made using eqn. (4.4). 

A A 

Smoothed estimates of ^11 ^‘^1^ selected 

frequencies 1 ■ 1, , , . , M, where M « N/2. This is to improve 
the reliability of estimates of Q(u)^). 

A 

Finally, is estimated for the selected frequencies co^, 

i*l, . . . ,M. 

The best procedure found for step 4 was one formerly applied 
to magnetotellur Ic studies, known as a "constant Q" analysis (Thayer, 
1975). After selection of "centered frequencies" u)^, 1 ■> 1, . . . , M, 

A A A 

M < N/2, the raw estimates of C.„(m ), C, ,(o} ) and ) are smoothed 

iz n 11 n zz n 

around each frequency o)^ using a Gaussian window W(a)^-u)^) 

W(«^-o)^) - (2/ir/a)^S) exp (--(w^-a)^)^/a!^^S^)) (4.7) . 

where S is the "selectivity." The selectivity defines the width of 
the Gaussian windows, which was set at S <■ 0.2 in this analysis. A 
"constant Q analysis" is one where the selectivity is held constant 
for all frequencies co^. This kind of window smooths less at longer 
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periods than shorter periods, which is necessary in order to v-taintaln 
relatively the same resolution at depth as near the surface of the 
earth. Specific listings of the programs used for this window are 
given in Thayer (1975). 

As noted earlier, smoothing Increases the reliability of the 
estimate. Statistically, it Increases the "degrees of freedom" of 
the estimate at frequency By increasing the window width, more 
data points are included in the smoothing and higher degrees of 
freedom are provided for the estimate. Since the constant Q analysis 
varies the window width with frequency, the number of degrees of 
freedom also varies with frequency. 

Within the "95% confidence interval" around an estimate of 

/N 

Q(o)^), there is a 95% probability that the true value of Q(u^) will 
exist. This confidence Interval depends upon the number of degrees 
of freedom. 

A 

For complex number Q(u)^) , the confidence interval is defined 

A 

in terms of a radius and a phase range Ai|)(u)^) where 


|q(o)^) I - r(to^) <_ IQ(w^) 1 1 IQ(w^) I + r(Wj^) 


4,(w^) - <_ <J)(w^) £ + A$(b)^) 


(4.8) 


where denotes "estimate." These are defined by Bendat and Piersol 
(1971). 


~-(h) 


°22<“l> 


A(^(o)^) ■ sin 


-1 


" r(u^) 


(4.9) 
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wh«r« 


n - number of degreee of freedom 

^2 n-2ta " percentage point of an F diatrlbutlon 
with n^^ ■ 2 and n2 ■ n-2 degreea of freedom 

/s 

^12^^!^ * **^^^^* coherency function betmen E((o^) 

and l(u^) defined by 


Y 




|Cij(Vl 




1 


An example of |q| and Q eatlmatea for March 5-9 » 1970 magnetic 
storm la preaented In fig. 4.1 a»b. 

To obtain a better estimate of Q(u>^)» the power spectra of 

/N A 

^12^^!^ and several magnetic disturbances were "stacked." 

Stacking accomplished by forming a weighted sum of spectra at 
each frequency o)^ from eight data sets listed in Table 3.1. Each 

A A 


spectral estimate of ^11 ^^1^ frequency for data 

set j was weighted with the factor n.(u3.)/In .(oi.) where n. Is Che 

jl jji j 

number of degrees of freedom for Chat estimate at frequency u) 1°^ 
data set j. The results for magnitude and phase are in fig. 4.2a 
and fig. 4.2b. 


To improve Che reliability of estimates of Q, power spectra 

2 

estimates at frequencies o)^ which had squared coherencies Y ^^2^^^!^ 
less than 0.60 were eliminated from the stacking procedure. This 


resulted in estimates of Q which have better reliabilities as shown 


In fig. 4.3a and fig. 4.3b. 
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Figure 4«1: Response function estimates for the March S, 1970 
magnetic storm. 

(a) Magnitude of Q 
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Figure 4.1: Reiponse function estlnatas for the >iarch 5, 1970 
magnetic storm. 

(b) Phase of Q 


i V 
f ' 
^ - 
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Figure 4.2: Response function for stacked data listed in Table 3.1 

A 

(a) Magnitude of Q 
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Figure 4.2: Response function for stacked data listed in Table 3.1 

A 

(b) Phase of Q 
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Figure 4.3: Response function for stacked data of Table 3.1 with 



/V 


(a) Magnitude of Q 



I 
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Figure 4.3: Response function for stacked data of Table 3.1 with 


2 

12 


> 


. 6 . 


(b) Phase of Q 
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A possible source of error In these estimates of Q comes from 

contamination of the total field measurements by the crustal anomaly 

field. Using a field model developed by Mayhew (1980), the crustal 

anomaly field is subtracted from the total field before separation 

into e(t ) and i(t ) contributions (Table 3.1). In larger magnetic 
n n 

disturbances such as of March 5, 1970, there is little change in 

e(t ) and i(t ) or in the response function. In smaller magnetic 
n n 

disturbances such as Sept. 29, 1969 there is a slight difference in 
estimates of Q(w) between the anomaly corrected data and the uncor- 
rected data (fig. 4.4 and 4,5). There appears to be no significant 
difference between stacked corrected and stacked uncorrected data sets. 
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Figure 4,4; Response funetion for September 29, 1969 magnetic storm 
with corrections for crustal magnetic aaomalles. 

(a) Magnitude of Q 
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Figure 4.4: Response function for September 29, 1969 magnetic storm 
with corrections for crustal magnetic anomalies. 

(b) Phase of Q 
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Figure 4.5: Response function for September 29, 1969 magnetic storm 
without anomaly corrections. 

(a) Magnitude of Q 




Figure 4.5: Response function for September 29, 1969 magnetic storm 


without anomaly corrections, 
(b) Phase of Q 
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CHAPTER V: INTERPRETATION OF THE ESTIMATED RESPONSE FUNCTION Q. 


In order to interpret the response function Q(ai) obtained from 

the satellite data (fig. 4.3), one must solve the "forward problem": 

given a conductivity profile, what is the response function? 

The response function Q is calculated for various conductivity 

profiles and compared to our best estimates of Q. All response 

4 6 

functions are calculated for periods 10 to 10 seconds which includes 
the range of our data. Chosen conductivity profiles include homogen- 
eous sphere, profiles by Banks (1972), Parker (1970) and other profiles 
with various upper mantle and near surface conductivities. 

V.l. Solution of the Forward Problem . 

The method presented here is a finite difference method operat- 
ing with real notation. 

Consider a sphere with some known conductivity profile depen- 
dent upon radius only. Surrounding the sphere is a nonconducting 
space with some source of time varying magnetic field. From Maxwell's 
equations, in MKS units 


V X H = 


§R. 

St 


+ J, 


(5.1) 


where magnetic field B * pH, D is the displacement field D - el and 
is the free current density. For slowly varying magnetic field, 
displacement current is negligible, inside of the sphere 


V X H * J. 

f 


(5.2) 
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Taking the curl of (5.2) anl using B * uH results in 


V^B » - pV X Jj 

Knowing aE * and V x E » - -g— , this becomes 

„2- 6B „ “ 

V B » ~ ^ E 


(5.3) 


(5.4) 


For constant conductivity, the 7o x E term will vanish. For conductivity 
that varies only in the radial direction, the B field produced by 
Vo X E term would only have 0 and (p components. Assuming that any 
initial charge density would decay to zero, Vo x E * 0 (Rikitake, 1966). 
So, the final equation to consider is 


2 - 

V^B 


pa 


6B 

6t 


(5.5) 


A complete solution to (5.5) is composed of a toroidal field 
T which would be restricted to within the sphere, a poloidal field S 
which would be detected outside of the sphere, and an arbitrary scalar 
field W which can be assumed to be zero. Since the eventual goal is 
to determine the field measured outside of the sphere due to both 
(5.5) and (5.6), the only solution of interest is the poloidal field S: 


S * - V X V X (rp) (5.6) 

So field equation (5.5) can be reduced to 

V^p = pa^ (5.7) 

6t 


where V x (rp) is the poloidal potential. 
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Outside of the conducting sphere, where no free currents exist 

V X H - 0 (5.8) 

Since (5.8) holds, B can be represented by a scalar potential U such 
that 

B - - VU (5.9) 

Since B of eqn. (5.9) consists of field contributions both of 
internal and external origin (relative to the sphere of radius a) , U 
can be represented by 

U - la (e^(t) (f)" + i^(t) (f)"""^^) (cos 0) (5.10) 

where a is the radius of the sphere, r is the variable radius, 0 is the 

magnetic colatitude, e is the external coefficient, i is the internal 

n n 

coefficient, and (cos 0 ) are the Legendre polynomials of degree n. 

n 

Assuming that the dipole term is sufficient to represent the time- 
varying disturbing field used in this study, U can be simply represented 
by 

U - a (e^(t) (|) + i^(t) (f)^pj (cos 0) (5.11) 

Since the response function Q(w) is 
F(i,(t)) 

Q(oj) - 

where F(ij^(t)) and F(ej^(t)) are the Fourier 
6j^(t), it is convenient to represent e^ and 


ideally 

(5.12) 

transforms of ij^(t) and 
i^ in terms of a Fourier 


series 
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“l ■ Flu““ "n‘ + “n‘ 

(5.13) 

1^ * Z Ii cos 0 ) t + I„ sin 0 ) t 
1 n In n 2n n 


wh.r. IB th. oxternal field contribution and 1^ Is the Induced field 

contribution to the potential. E, , E„ , I, , I_ are the Fourier 

in zn in zn 

coefficients to the power series. Substituting (5.13) into (5.11) gives 
a simple power series representation to the potential U, 

From the poloidal field (5.7), p can also be represented by a 
Fourier power series using the dipole term (cos 9 ) 


p “ r ^EA„(r) p!^(cos 0) cos 0)_t + C„(r) P^(cos 6) sin w^t (5.14) 


n n 


n 


Substituting (5.14) into the differential equation (5.7) results in a 
coupled pair of equations 


d^A 


dr 


" O -2 A 

■» - 2 r A 
^ n 


- w C 
p n n 


,2 

d C 


dr 


I -2 r“2 c 

^ n 


1 

— 0) A 

p n n 


(5.15) 


where P * is resistivity and y * is the free space permeability. 

(We assume y * y^ for the sphere.) 

Using the conditions at the surface of the sphere (r = a) for 
continuity of the B normal and H tangential components, we arrive at 
the boundary conditions at r = a 
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»3rl. + 

a _E I + C„(n) •- E_ 

dr ' a " 


(5.16) 


Since these equations (5.15) to (5.16) are solved numerically, 
It Is desirable to express them in nondlmensional form; we make the 
following change of variables 


r « a - (1 - z) N(S 
2 

- a K 

n n 

2 

C - a C 
n n 


(5.17) 


Substituting (5.17) into (5.15) yields (5.18) where c ■ N6/a and 
6 - t2/yw0)^^^ 


2^ 

d A 2 2 ^ 2 - 

2_c ^ ^ . 2 C 

dz^ (c(z-l)+l)^ " " 


2^ 

d C ^ 2 ^ 

n 2 n 2 

— ^ £_£ c -- 2 A 

2 n n 


dz (c(z-l)+l)‘ 


(5.18) 


with boundary conditions at z ■ 1 


A 



+ c A (1) 
n 




A 

dC 
n 

dz 


z»l 


+ c C (1) » c E 
n 2 


2n 


(5.19) 
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where c ~ . 
a 

Examining the pololdal field, 5 must go to zero at r * 0, so 

A (z*0) and C (z ■ 0) goes to some small number asz-^0, since the 
n n 

strength of the inducing field decreases with depth. By normalizing 
with "N*' skin depths, z»0 implies that the field is very small when 
it has penetrated N skin depths. So, additional boundary condltons at 
z * 0 are 


A^(0) - a 

(5.20) 

C^(0) - 6 


where a and 3 are small and are assigned according to "N". 

The actual solution of (5.18) is accomplished using a Runge' 
Kutta finite difference algorithm, for 4 simultaneous first order 
differential equations 


Yi - A^(z) 

A 

Y ■ A ' 

^2 \iz) 

A 

Y„ - C (z) 
■i n 


then (5.18) becomes 
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dY, 

dz ^2 


2c ^ 


F„ * 


(c(z-l)+l) 

dY» 

— i B Y 
dz 


1 2 

+ 2 Y3 


dY^ 

? m 2 

4 dz 


2ch 


(c(z-l)+l) 


^ - 2 Vj 


(5.21) 


where c ■ and z is radirl distance, 
a 

A theorem for differential equations (Boyce and DiPriroa, .1969) 
states that for a system of n first order, linear, homogeneous differ- 
ential <■ .vations 


P^j(t)Xj + 




each solution x ■ 4>(t) of this system can be expressed as a linear 
combination of x^^\ ... , x^*^\ with arbitrary coefficients c^ 


X » c-x^^^(t) + ... + c x^”^(t) 
1 n 


in exactly one way. Knowing this, 4 independent solutions to (5.20) 
are obtained using arbitrary A^(0), k^(0) , C^(0) , C^(0)> The final 
solution, (z) + + C3 A„,(z) + c,A„,(z), must meet the 

four boundary conditions (5.19) and (5.20). This requires solution 
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to a set of four equations with four unknowns Cj , c^. c^, c^j 




^ ‘=2^na^°> ‘^3^n3<°> ‘=4‘^n4<0> - 6 


c(3/2 - Cj^A^j(l) ~ C2A^2^1) ■" *^3^03^^^ “ Vn4^^^^ 

‘=<3/2 E2„ - c^C^^(l) - C2C„2<1> ' Vn3<^> “ ‘=4^n4<^» 


(5.22) 


N6 


fl t If 


where ^ and £ 2 ^^ are Input amplitude (eqn. 5.13) and 

denotes first derivative with respect to z. 

A A 

Once the final solution for An(z) and Cn(z) is obtained, for 
particular conductivity profile and frequency w, the response function 
for the assumed and £ 2 ^^ input (arbitrary since Q Involves a ratio, 


so » E 2 j^ ■ 1) is computed 


«<V 


®ln - 3 ^2„ 


(5.23) 


where j ■ /-I and and 12 ^ are derived from the boundary conditions 
(5.19). The response magnitude |q| and the phase (p(Q) can be 
expressed by 


iQl 


<^L 

«ln + 


‘V2n-Vln> 

(})(Q) - arctan TiT-eT-TTTEo ) 
^ In In 2n 2n-^ 


and 


(5.24) 
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Alt«rtuitive methods for solving the forward problem Involve 
solving the "R" part of the pololdal field equation (see eqn, 5,7) 
after separation of variables. One method requires a series of shells 
of constant conductivity where the analytical solution for each shell 
Is known (Banks, 1969). The Runge-Kutta method can also be applied 
to the "R" equation In complex notation (Parker, 1970), Eckhardt 
(1963) transforms this same "R" equation Into a nonlinear first order 
differential equation and solves It numerically. 

V.2. Behavior of the Forward Problem Solution . 

A 5 

Response functions were calculated for periods 10 to 10 
seconds, mantle conductivities .001 to 100 mho/m and an accepted 

5 

core conductivity of 3 x 10 mho/m (Stacey, 1977). This highly 
conducting core also Insured that the boundary conditions at N skin 
depths would remain within the dimensions of the earth. 

The response function for the range of periods considered is 
most sensitive to the upper regions of the earth. Further, the degree 
of sensitivity depends upon the conductivities Involved and the rate 
of c^uage of conductivity with depth. Large conductivities tend to 
hide variations in conductivity while lower conductivities are more 
transparent. Large changes in conductivity over small changes in 
depth are also more visible than gradual changes. For example, Parker's 
1970 profile (fig, 1.2) compared to a homogeneous sphere of .1 mho/m 
shows only a small change in response function (fig. 5.1). Parker's 
profile is a smoothly varying model with the upper mantle having a 
conductivity of .1 mho/m. In con-trast. Banks 1972 model (fig. 1.1) Is 
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Figure 5.1s Parker's 1970 model compared to homogeneous sphere of 
.1 mho/m, 

(a) Magnitude of Q 



Figure 5.1: Parker's 1970 model compared to homogeneous sphere of 
.1 mho/m. 

(b) Phase of Q 


PHASE OF Q VS. PERIOD 



LOG PERIOD 
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composed of a series of steps with large changes in conductivity and 
an upper mantle conductivity of .01 mho/m. Banks' model compared to 
a homogeneous sphere of .01 mho/m (fig. 5.2) shows a large contrast 
in response functions. 

In assigning starting depths, it was found that the calculated 
response function was fairly insensitive to changes in N (number of 
skin depths) for N > 5 and a and 0 of order .001. This is particularly 
true for conductivity profiles with near surface values greater than 
.01 mho/m. Banks (1972) also calculated the response function for his 
conductivity profile and for Parker's (1970) using the analytical 
solutions for constant conductivity shells. Our response functions 
agree fairly well with Banks' calculation of Parker's profile, but 
disagree with Banks' calculation of phase for Banks' profile. Our 
phases for small conductivities tend to Increase more rapidly with 
period than do Banks'. Low conductivities have large skin depths, 6, 
so "c" (5.17) is greater than one. The phase is greatly affected 
by how much of the layer with c > 1 is included in the finite 
difference calculation. The phases presented here are the largest 
allowable. However, the magnitude calculation is essentially exact. 

General behavior of the response functions is best illustrated 
by calculating Q for homogeneous spheres (fig. 5.3a,b). As conduc- 
tivity increases, the magnitude of the response function, |q|, 
increases while the phase, (j), decreases as a function of period. 


I 



Figure 5.2: Banks’ 1972 model compared to homogeneous sphere of 
.01 mho/m. 

(a) Magnitude of Q 




Figure 5.2: Banks’ 1972 model compared to homogeneous sphere of 


.O'v mho/m. 

(b) Phase of Q 


PHASE OF 0 VS. PERIOD 



LOG PERIOD 





68 


Figure 5.3: Response functions of homogeneous sphere of varying 
conductivity. 

(a) Magnitude of Q 


MAGNITUDE OF 






Igure 5.3: Response functions of homogeneous sphere of varying 
conductivity. 

(b) Phase of Q 


vs. PERI 
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V.3. Comparing Eatlmatad Q to Calculated Q, 

As a first stop in interpreting the estimated values of Q> 

calculated response functions for simple conductivity profiles were 

compared to the estimated response. When response functions of 

homogeneous spheres are compared to estimated responses (fig, 5.4), 

the magnitudes lie in the vicinity of the .01 and .02 mho/m curves 

5 

for periods greater than 10 seconds. For periods less than 10 
seconds, the estimated magnitudes lie well below the .01 curve. For 
all periods considered, the phase lies below or within the 95% confi- 
dence Interval of the .1 mho/m curve. Neither estimated magnitudes 
nor phase show the same mono tonic shape as do the response function 
curves of homogeneous spherec. Clearly, a homogeneous sphere of any 
conductivity cannot explain this data. 

Next, the estimated response function was compared to the 
response function for Banks' 1972 profile (fig. 1.1, fig. 5.5a) and 
Parker's 1970 profile (fig. 1,2, fig. 5.5b). The estimated magnitude 

(fig. 5.5c) is compatible wij:h that of Banks' for periods greater 
5 

than 10 seconds, but the phase is not (fig. 5.5d). The phase curve 
for Parker's model lies within the 95% confidence interval of the 
estimated phases but does not show the same monotonic shape. In all 
of the foregoing comparisons, the phase and magnitude estimates seem 
to indicate conductivities which are nearly an order of magnitude 
different from each other. Since the phase estimates contain non- 
physical negative values, more credence is given to the magnitude 
estimates. Modifications of Banks' conductivity profile (fig. 5.5a) 


Figure 5.4: Comparison of best estimate of response function to 


those of homogeneous spheres , 
(a) Magnitude of Q 
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Figure 5.4: Comparison of best estimate of response function to 
those of homogeneous spheres. 

(b) Phase of Q 




Figure 5.5(a): Banks' 1972 conductivity profile as used in this 


study. 
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Figure 5.5(b): Parker's 1970 conductivity profile as used in this 


study . 
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Figure 5.5(c): Comparison of magnitude estimates to those of Banks’ 


and Parker. 
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Figure 5.5(d); Comparison of phase estimates to those of Banks and 


Parker . 
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should yield response functions that are compatible with estimated 
magnitudes. 

Since the frequency range considered should be sensitive to 
near-surface conductivity structures, layers of varying conductivity 
were superimposed onto the top 500 km of Banks' model. With a 30 km, 
.001 mho/m layer replacing the top 30 km of Banks' model, the calcu- 
lated response magnitude agreed well with the estimated magnitude 
for periods less than 10^ seconds (fig. 5.6, fig, 5.7a). Changing 
Banks' .01 mho/m layer to .02 mho/m, but maintaining the 500 km lower 

boundary (fig. 5.7) in this model, improved the agreement of esti- 

5 

mated magnitude to the calculated magnitude in the vicinity of 10 
seconds; the values for longer periods were too high. This establishes 
a conductivity range of .01 to .02 mho/m for depths just greater than 
30 km (fig. 5.7). To test the sensitivity to layer thickness, the 
lower boundary of the .01 layer was varied from 500 km to 350 km 
(fig. 5.8); the magnitude shows very little change in this frequency 
range. Variance of conductivity from .01 to .02 mho/m with a lower 
boundary of 350 km (fig. 5.9) also shows very little change. 

Changing the conductivity from 2. mho/m to i.mho/m for the layer below 
500 km did not change the response function at all. 

V.4. Ocean Effects on Response Function Q. 

An electric current travelling a near-surface e.quitorial 
path would encounter continents with a lower conductivity than the 
oceans. The effective conductivity of an ocean-continental medium 
can be considered as a simple model in which the oceans and continents 
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ar« slabs of uniform cross-section but have a variable thickness. 

From Ohm's law one can calculate the effective conductivity 
of a medium made up of slabs of varying thickness and two conductlv- 
ltl« Oj and a^: 


a 


h ^2 


j 


T. t 


11 


"2 


(5.25) 


th 

where t^^ is the thickness of the 1 slab of conductivity and 

til 

is the thickness of the j " slab of conductivity 02* T is the 

total thickness and all slabs have the same cross-section area. 

Applying equation 5.25 to an ocean-continental medium, 

surface area fractions of ocean t^d/T) and continent (I t„./T) 

i 11 j 

were estimated along the equator between 30" N and 30" S latitude. 

Using a 5“ X 5" grid, estimates were made for each 30" longitude 

spread along the equator. Using typical continental conductivities 
-3 -1 

of 10 to 10 mhos/m derived from magnetotelluric data (Gough, 
197A) , effective oceanic conductivities were calculated. For an 
ocean of 3.3 mho/m, these values are: 


-3 -3 

3.78 X 10 mhos/m for 10 crust 

-2 -2 
3.76 X 10 rahos/m for 10 crust 

3.49 X 10 ^ mhos/m for 10 ^ crust 

The above relationships indicate that the effective conductivity is 
dominated by the continental conductivity. 

Response f'mctions were calculated for effective ocean 


Figure 5.6; Comparison of estimated magnitudes to calculated magnitudes 
for Banks’ 1972 profile and for Banks’ profile with a 30 km 
surface layer at .001 mho/m. 




Figure 5.7: Effects of variation of second layer (lower boundary at 
500 km) from .01 mho/ra to .02 mho/m. 

(a) with .01 mho/m layer 

(b) with .02 mho/m layer. 
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Figure 5.7: Effects of variation of second layer (lower boundary at 
500 km) from .01 mho/m to .02 mho/m. 

(c) Response magnitude estimates compared to calculated 
response magnitudes. 
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Figure 5,8s Effects of variation of second layer thickness. 

(a) Lower boundary of .01 mho/m layer at 500 km. 

(b) Lower boundary of ,01 mho/rn layer at 350 km. 
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Figure 5. 


8: Effects of variation of second layer thickness. 

(c) Response magnitude estimates compared to calculated 
response magnitudes. 




MAGNITUDE OF Q VS PERIOD 
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Figure 5.9: Effects of variation of second layer (lower boundary at 
350 km) from ,01 mho/m to .02 raho/m. 
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conductivitiss superimposed onto the model (fig. 5.7). The 
.00378 mho/m ocesn of 3 km thickness made virtually no difference 
In calculated response function when compared to the original model . 

The .0376 mho/m ocean raised the response function magnitude for 
periods Just below the data range. For effective ocean conductivities 
greater than .0376 mho/m, the response function magnitude is signifi- 
cantly raised in the frequency range of Interest. Since our estimated 
magnitudes are much lower than these values, an upper bound of .01 mho/m 
is placed on continental conductivity for the top 3 km (fig. 5.10). 
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Figure 5, 10 5 Conductivity profile for estimated response function. 
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CHAPTER Vis TEMPERATURE IN THE MANTLE 

VI. 1. Relation of temperature to conductivity 

Electrical conductivity estlinates for the earth are often 

related to temperature. This Is a valid relation if the conductivity 

is Intrinslcly controlled by temperature and if the composition is 

known. Estimated electrical conductivity for the upper mantle and the 

-9 2 

crust is within the range of semiconductors (10 to 10 mho/m). For 

these materials, conductivity depends upon populating the "conduction 

band" (fig. 6.1) with electrons from the "valence band". In the 

intrinsic temperature range, electrons are thermally activated into the 

conduction band, where they become mobile and contribute to the 

conductivity. This temperature dependence is expressed by 

0 ■ a exp (-E/2kT) (6.1) 

o g 

where k is Boltzman's constant, T is temperature, 0^ is a constant 

characteristic of the semiconductor, and E is the activation energy 

S 

or "gap energy". A plot of the logarithm of resistivity vs 1/T 
i.*eveals a straight line (fig. 6.2) if the semiconductor is intrinslcly 
controlled. 

At higher temperatures, a semiconductor may become an 
"ionic conductor". Hera, the charge carriers are ions (lattice defects) 
rather than electrons (electron-hole pairs) . Conductivity is controlled 
by the rate of diffusion of the ions. For sufficiently high temperatures, 
ion diffusion is also temperature controlled, implying that conductivity 
can be represented by 


¥ 
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Figure 6.1: In an intrinsic semiconductor, electrons are thermally 

excited from the valence band into the conduction band, where they 
become mobile. (Kittel, 1971) 
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Figure 6.2; A material may be aa extrinsic semiconductor at low 
temperatures and an intrinsic semiconductor at higher temperatures. 
Here, germanium is doped with varying amounts of galliurfl resulting 
in different temperatures for intrinsic semiconduction to occur. 
(Kittel,1971) 
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a - 0^ exp (-E/kT) - Nqy exp (-E/kT) (6.2) 

where “ Nqy, N Is ion concentration, q is charge, y is ion mobility, 
and E is the activation energy for an ion to move from one location 
to a similar one. 

Pressure effects probably have little or no Influence 
on intrinsic semiconduction where the charge carriers are electrons. 
Experimental results for pressures up to 8 kb indicate that conductivity 
is weakly affected (Duba,et.al. ,1974) . Pressure effects might become 
more significant where ionic conduction is the dominant mechanism 
since lattice compaction may affect ion mobility. However, it is 
difficult to say what such pressure efects might be (Misener, 1973) . 

In any case, it is probably safe to assume that temperature effects 
will dominate over pressure effects. 

VI. 2 Composition of the mantle 

The velocity ot the seismic P-wave in the upper mantle is 

3 

8.2 + 0.2 km/ sec which indicates a density of 3.3 g/cm . This 
observation along with petrological considerations, suggests an upper 
mantle composition of peridotite (olivine + pyroxene + spinel) and/or 
eclogite (pyroxene + garnet + quartz). Eclogltes, which are thought to 
be of mantle origin, tend to have slightly higher density than expected 
for the mantle. Peridotite is generally considered the more likely 
mantle material. 

A peridotite mantle is also suggested by the common 
occurrence of peridotite inclusions found in Kimberlite pipes. 

Among these xenoliths, olivine principal mineral, 
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with an averaga ollvlne/orthopyroxene ratio of 2/1 (Ringvrood, 1975). 
Garnet Iherzolite — the moat conunon peridotitic inclusion — is com- 
posed of olivine 64%, orthopyroxene 27%, clinopyroxene 3%, pyrope - rich 
garnet 6% (Ringwood, 1975). Such petrological observations lead to the 
assumption that the properties of the upper mantle can be approximated 
by the properties of olivine This single mineral model Js 

further supported by Birch (1969), who concluded that the elastic 
properties of the upper 2C0 km of the mantle could be represented 
to first order approximation by a homogeneous layer of olivine. 

A peridotitic mantle, composed mainly of olivine, is also 
a convenient model for explaining various selsmologlcal occurrences. 
First, the almost worldwide increase in seismic velocity at the 
Mohorovlcic boundary could be attributed to the transition from a 
gabbroic lower crust to a peridotitic upper mantle. Secondly, the 
increment in seismic velocity observed at 420 km could be explained 
by a phase transition of olivine from a low-pressure orthorombic form 
to a high pressure spinel structure. Finally, the low-velocity zone, 
occurring between 60 and 150 km beneath the oceans and in the vicinity 
of 200 km under shield areas, is believed to be due to Incipient 
partial melting. This localized phenomenon could be the result of the 
depressed melting temperature of peridotlte at 20-30 kb in a CO 2 and 
H 2 O environment (Eggler and Kushlro, 1979; Fig. 6.3). 

Alternative mantle models have Included an eclogite mantle 
and an eclogite + peridotlte mantle. As already mentioned, an eclogite 
mantle would have a higher density than would be expected from seismic 
observations. Also, an eclogitic mantle would require an intermediate 
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Figure 6.3; P-T diagram for the reaction En + Mag • Fo + CO^. 

(A) Calculated boundary for the reaction. (B) Shovm are the phase 
assemblages for a peridotite composition in the system Mg 0 -Si 02 “C 02 “ 
H.^0 containing less than 20 weight percent 002* The solidus is for 
peridotite with amounts of K 2 O and CO^ small enough to be buffered 
within the system. (Eggler and Kushiro,1979) 
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Mg0-Si02-CCD2-H20 
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step between the gabbrolc crust at 15 kb and the ocloglte mantle at 
pressures above 20 kb (Carmichael, Turner, Verhoogen, 1974), A gamet- 
granullte assemblage (garnet 4 pyroxene + plagloclase) could exist 
at the 20 + 3 kb pressure regime but, the complex "phase transition" 
of gabbro-granulite-eclogite is incompatible with the sharp toho 
discontinuity. Although an eclogite mantle is favored by some for its 
compatibility with certain magma genesis models, a peridotltlc mantle 
seems the preferred model. 

VI. 3 Electrical conductivity of olivine ^FOg 8 _ 94 ) 

Early investigations (Mlsener, 1973) Indicated that olivine 

(FOgg^g^) is an extrinsic (impurity controlled) semiconductor below 

800°C. Later measurements (Duba,1974) showed the expected temperature 

dependence for temperatures greater than 800^C. Because of kinetic 

effects and apparatus problems, accurate measurements of conductivity 

is limited to temperatures below 1500°C. 

A necessary condition in all high temperature experiments 

is to have controlled oxygen fugaclty (f^ ). Nltsf«n (1974) showed 

"2 

that olivine (herein chosen to be the representative case of 

is stable over a limited range of oxygen fugacitles. Outside of 

the stability field, olivine wi.ll oxidize or reduce to other minerals. 

Duba (1976) found that electrical conductivity varies less than 1/3 

of a magnitude for changes of f_ within the stability field. However, 

^2 

for f^ variations outside of the stability field, conductivity can 
^2 

change by several orders of magnitude (fig, 6. 4). These reaction 
products apparently provide better electrical pathways. 
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Figure 6. A; Stability field of olivine (Fo^q) defined with 
regard to temperature and pressure (Nitsan» 1974) . Also included 
are buffers OQPM (fayalite-quartz-pyroxene-magnetlte) , OQM (fayalite- 
quartz-roagnetite) , MW (magnetite-wustite) , and WI (magnetite-iron). 
"Error bars" are the region where Duba (1976) tested the variation 
of conductivity with variation in f_ within the stability field and 
found that conductivity varied less than 1/3 magnitude. 
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It ii r«a«onabl® to «xp«ct r.h« upptr mantl* to b« in 
tharmodynamlc •quillbrlum. It X» unlikaly that laboratory axparimants 
totally raach aquilibrium bacauaa of tha briaf tima that thay axict. 

This idaa la conflrmad by Duba'a (,at.al. , 1974) obaarvation at flxad 
tamperatura, praaaura, and whara conductivity variad with tima 

U2 

ovar a period of a few houra. 

Conaidarlng all of thaae axparimental dlfflcultiea, 
conductivity for two aamplas of olivine ('^Fo^q) have bean maaaurad and 
are praaauted in figure 6.5 (Duba and Nicholla, 1973j Duba, at. al., 1974, 
Duba,1976). For coraparf.aon, the conductivity of baaalt (Duba,et.al. , 
1975) is Included (fig. 6.6), but these results should be viewed with 
caution due to the inherent experimental difficulties encountered when 
dealing with a rock rather than a single crystal. 

As previously mentioned, pressure effects on the Intrinsic 
conductivity of olivine is probably small for depths less than 400 km. 
However, pressure could change the shape of the stability field with 


regards to . This change can be estimated by (Carmichael, et.al. ,1974) 

log(f ) - lug(f ), . _ - AV^^(P-l ) (6.3) 

^2 ^ ^ 2.303 RT . 

where P ia pressure, R is tu., gas constant, T is temperature, and 


is the change in volume for the solids in the reaction. For the reaction 




Figure 6.5: Variation of conductivity as a function of temperature 

for single crystal olivine ('VFo^q) . (a) Olivine from San Carlos Indian 
Reservation (Duba and Nicholls,1973;Duba,1976) (b) Olivine from St. 
John's Island (Duba, et.al.,1974; Duba, 1976) All are under controlled 
oxygen fugacity. 




Figure 6.6: Flectrical conductivity as a function of temperature 

for basalt at controlled oxygen fugacity (Duba,et.al. , 1975) . 
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become log(f ) « -7.94 at 32 kb and log(f- ) ■ -6.02 at 124 kb. The 
"2 ^2 

rate of change of log(f_ ) with pressure is different for different 

^2 

reactions, so the overall shape and location of the stability field 

will change with pressure. Compared to the stability field at 1 bar, 

there will be little change at shallow depths (100 km, 32 kb) and some 

change before the transition zone (400 km, 124 kb). Still, since the 

electrical conductivity does not change drastlcly for within the 

2 

stability field, this pressure effect Is of minor Importance In this 
study. 


VI . 4 Temperature as a function of depth 

Numerous chemical analysis of basalts from around the world 

Imply that . the f_ within the mantle should be near the fayallte- 

magnetite-quartz buffer (fig. 6. 7). The olivine conductivities (fig. 6. 5) 

used here are at controlled f_ ’s within the stability field of 

“2 


olivine (FOgg) but vary from the oxidizing side to the reducing side 

as a function of temperature (Duba,1976). Knowing that variations In 

f_ within the stability field result In changes of conductivity up 
^2 

to 1/3 magnitude, a certain amount of error Is allowed for f^ . Also, 
knowing that the depths obtained In the estimated conductivity profile 
(fig. 5. 10) derived from "Q" are upper limits, error is allowed for 
location. Considering these error contributions, a temperature profile 
based on the conductivity of olivine <''^<Fo^q) and basalt Is presented 
(fig. 6. 8). The two profiles differ a great deal, but this difference 
may largely be due to experimental errors in the determination of 
conductivity in basalts. . 
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Figure 6.7: Oxygen 

the world. There is a 


fugaclty calculations for basalts from around 
strong tendency to follow the fayallte-magnetlte' 


quartz (FMQ) buffer (Haggerty, 1976). 
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Figure 6.8: Temperature profiles estimated from response function 

Q and associated conductivity. Conductivity- temperature relation for 
olivine is from fig. 6. 6. Conductivity-temperature relation for basalt 
Is from fig. 6. 7. 
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Teapcratura •■tlnatci for depths down to 200 km were made 
using detailed compositional analysis of kimberlites (Boyd, al. , 
1973; fig. 6.9). The resulting temperature profiles (fig. 6.10) inter- 
sect ths ”Q” profile (based on electrical conductivity of olivine) in 
the 150 to 200 km range. For shallower depths, the geochemical 
geotherms are much lower than what is indicated by Q. This disagreement 
may be due to several factors. One very likely error is in assuming that 
olivine is the controlling agent for shallow conduction. Other materials 
are probably controlling conductivity at shallow depths. At temperatures 
below 800°C, conductivity is impurity controlled rather than temperature 
controlled — leading to conductivity values higher than would be 
expected. Finally, very low conducting surface layers (10”^mho/m or 
less) would not be seen by the Induction method employed here for 
frequencies of 0.2 to 2.0 cycles/day (according to tests employing 
the forward problem calculation for Q). The pyroxene geotherm combined 
with the olivine measurements indicate a conductivity of lO'^mho/m 
near the surface — this value is too low to be detected. 
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Figure 6.9: Geotherm estimates based on detailed compositional 

analysis of Iherzolite nodules from kimberlites. Slight differences 
In geotherm will occur depending on how pressure is calculated 
(Boyd and Nixon, 1973). 




Figure 6.10: Comparison of **Q" geotherm to pyroxene geotherm of Boyd 

and Nixon (fig. 6, 9). Q geotherm Is based on electrical conductivity 
of olivine (fig. 6. 5). 


I 'KW 








129 


Figure 6.11: Expected electrical conductivity profiles for the upper 

250 km if pyroxene geotherm is correct (fig, 6. 9). Conductivity 
estimates are based on olivine (fig. 6. 5) and basalt(fig. 6.6) . 
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CHAPTER VII: CONCLUSION 


A global electrical conductivity prol'lle for the upper mantle 
and crust Is obtainable through analysis of satellite magnetic field 
data. The resulting profile Is based on estimates of Q(w) for 
frequencies of 0.2 to 2.0 cycles/day. This range Is higher than that 
generally used In global Induction studies using land observatory 
data. 

The upper mantle has a characteristic conductivity of order 

_2 

10 mho/m. This result Is compatible with Banks’ (1972) model but Is 
not compatible with Parker's (1970) values of order lO^^mho/ra. 
Considering laboratory measurements of conductivity of olivine (Duba, 
1976), Parker's conductivity would correspond to temperatures that 
would melt the mantle. 

Shallow structures are Indicated by the response function 

^ —2 

Q(u). An upper limit of 10 mho/m Is placed on the top 3 km of crust. 

_3 

Beneath this, a conductivity of 10 mho/m extends downward to at 

—2 

least 30 km depth before Increasing to values of order 10 mho/m. The 
_2 

bottom of the 10 mho/m layer Is difficult to define because of the 
limited range of frequency, but the lower boundary has an upper 
limit of 350 km. 

In relating temperature to conductivity, the temperature 
profile derived from conductivity estimates Is acceptable for depths 
greater than 150 - 200 km. For shallower depths, temperatures based 
on olivine conductivity are too high for a perldotltlc mantle. 
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This disagreement between the "Q" geotherm and the pyroxene geotherm 
leads to several possible problems. One is in assuming that olivine is 
the contolllng agent in shallow conduction . We know from numerous 
geological studies that it is not the main material in very shallow 
depths. Basalt conductivity measurements may be more appropriately 
applied to this study, but laboratory measurements tend to be 
unreliable. Another problem concerns the sensitivity of this induction 
method to shallow, low conducting layers. The pyroxene geotherm combined 
with the olivine measurements would Indicate conductivities of order 

_5 

10 mho/m — such conductlvites would not be seen by Q estimates for 

frequencies 0.2 to 2.0 cycles/day. Whether or not the surface conductl- 
-4 -5 

vlty is 10 or 10 is a moot point with regards to temperature since 

at such low temperatures, the conductivity is Impurity controlled. 

As for the conductivity profile itself , magnetotellur ic 

measurements are compatible with the shallow structure — but local 

field measurements also' Indicate values both higher and lower than 

_2 

estimated from Q. However, the upper mantle value of 10 mho/m seems 
quite reliable. 
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Appendix A: Polar Orbit 


(Illustration on next page) 
Characteristics of satellite orbits; 


Satellite 

launch date 

Inclination 

perigee 

OGO 2 

Oct. 14, 1965 

87.3° 

410 km 

OGO 4 

July 28,1967 

86.0° 

410 km 

OGO 6 

June 5,1969 

82.0° 

400 km 


(Langel,1979) 


apogee 

1510km 
910 km 
1100km 



Appendix Bs Main program for finding e(t) and i(t) 
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c 

C 


400 

501 


MAME IS 'FIPIU* 

•• AS OF MAHCH 7.19B0 

ATTEMPTS TO WXITE Ow TEMPORAMV OISK FILE. 

WOTF .ICL 

CALLS fiUESSS FDK INTEHMOLATIOM 
CALLS ECOMF FOR A.MfJi-ALV COKHECTIOM 

PURPOSE IS TO calculate EITI.KT), 
nRIGIN'AL PROfiKA). WAS 'TS?' WHICH CALCULATES 
E(F».I|F|. RATIO. AMO PLOTTINC, ROUTINES. 

lMTBr,ER*2 INIIO.SOI.IPASS 
dimension AREA(3f>l i.MSIl'ISn) 
dimension EIMATnSO.SI 
dimension 0UT2(150.4) 

COMMON /FLOCOE/ST.CT.SPH.CPH.R.mmax.HT.HP.RH.R 
DATA RATIO/. 10/. RAD/57. 29577V/ 

DIMENSION P0T(5.15O) 

Ht-AL^R OTAPE 

DATA DTAPE/'0S420(S'/. SEAS/1. /.lOTVPF/l/.MlO'IT/H/. IPASS/ 0/ 

DATA IO/O/.NFILE/I/ .IX/0/ 

EOin VALENCE! IN! 1, 1 I .AREA! 52 » ) . !M.JD,AREA! 1 1 >,! MSU'! 1) .AREA! 2) ) 
COMMON /CAL. SWITCH 
BRHBEIMAT FOR E ANO I: 

SWITCH-0.0 

SWITCH-l.n INDICATES CALCOMP PLOTS 

ITEST-0 

IE-0 

L-1 

call PIELDC.CO. .0..0..1V(Sfl.,lfl.L. A.D.C.FI 
rewind 5 

READ ! 15.4001 NJDS.NJOE.MSS.MSE. IPR 
WRITEI7.400) NJDS.NJOE.MSS.MSE. IPK 
FORMAT (215.31101 

WRITE (B.SOn NJDS.NJOE.MSS.MSE .IPR 
FORMAT (2X.2I5.3I10) 

TM«CNVJUL!NJDSI 
LREAD ■ 1 
LPRINT-1 
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1 


6 


A 


C 

C 


NHX * 14 

CALL MOUNT! lOTVPEfNUNIT.OTAPE) 

INC-0 

OBM - 0. , 

IBS ■ 0 

INL-0 

IPC-0 • 

TT-0.0 

RA-0.0 

IC-0 

DUMMY VARIABLE 1C USED fOR PLOT SUBROUTINE. HAVE 
MOOIEIEO ORIGINAL PROGRAM BY ADDITIONS AND BY 
calling other SUBROUTINES AT NUMBER 500 INSTEAD OF STOPS. 

READ 05. 4501 TO . , ^ 

WRITE«7,450» TU 
FORMATOX.F5.2) 

WRITE(6,470I TO 

FORMAT! IX. 'LOCAL TIME USEOt • .2X.F5.2 I 

A-TU* l.O 

BA-TU-l.O 


call POSN! IOTYPE.NUNIT.NFILEI 
CALL FREAD!AREA.NUNIT.LEN.C500.C5I 
IF!LEN.NE. 12041 CALL ABEND!5I 
IF!MJD.LT.NJDSt 60 TO 5 
IF!MJD.6T.NJ0EI GO TO 500 
INC « INC ♦ 1 
PO 50 I- 1.50 

IFIMSIN! n.E0.-V99> GO TO 5 
IF!MSIN! n.LT.OI GO TO 50 
IF!MJO.GT.NJDS.ANO.MJO.LT.NJDEI 60 TO 8 
IF!MJO.EO.NJOS.ANO.MSIN!I).LT.MSSt 60 TO 50 
1F!MU0.E0.NJ0E.AN0.MSIN! n.GT.MSEI 60 TO 500 
IF! IPASS.EO.OI GO TO 9 
1F!IPASS.E0.1N!B«IM GO »0 
IF IIP.LT.20I 60 TO , - 

IF!DPM.GT.5. I GO TO 9 

IF!A.LT.24.0.AND.BA.GE.O.O» GO TO 12 < 

I HAVE REPLACED .OR. WITH .AND. IN THE ABOVE STATEMENTII 

LOCAL TIME CHOICE CORRECTED HERE. I HOPEI I III I 

IF!TML.6T. 0.0. ANO.TML.lt. 1.01 60 TO 2 

IF!A.6E.24.0) 60 TO 3 

BA- BA4 24.0 

RA-ABS!TML-BAI 

BA-BA-24.0 

IF!RA.6T.1.0) 60 TO 9 
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r,n TQ 15 

9 RA«ABS(A-TMi.) 

lF(HA.r,T.2.0) GO TO 
GO TD 15 

2 1PIA.GE.2A.O) GO TO 
RA«AHSITML-*«AI 
IF(‘<A.r,T.2.0l GO TO V 
GO TO 15 
1(S A-A-2A.0 

RA«ArtS« A-TMi. » 
timA*?A »0 

IP(KA.GT.l.O) GO TO V 
Git TU 15 

1? H‘(TMi.,i.T,RA,nR.T^L,GT.AI GO TO R 
15 Ff S««S 

UTHK • Rl>’S/360000«. 

ISw ■ 1 

n on ion <»1,IP 

lf=<Pr)T(J,K).LT.A5.n.DH.POT(l.K(.GT.n5.) GO TO 100 
1* Urn PnT(?,<» f <S371. 

Th > ( 90. - P0T(4.K » »/RflO 
SI - S1N(TH» 

CT ■ f,iiS(TH) 

PH ■ P0T(5.'< l/RnO 
SPh ■ .SI^'(PH) 

Oph ■ COS(P«) 

CALL PIELO 
A'M ■ A571./P 
AHR3 ■ AHKW*? 

THO « P0T( l.Kt/KAn 
STD » SIN(THD) 

CTO ■ CnS(THD» 

CnSA » (.9H027 - CT«»CTO 1 /STD/ST 
lP(CltSA.GT.l. ) COSA - 1. 

if(Chsa.lt,-i. ( COSA • -i. 

SI'9A ■ SORT (1. - CnSA**?» 

BTD • HTOCnSA - HP*Sl^fA 
RPD « MPOCOSA ♦ HT*SH^A 
IPdSM.EO.U GO TO 1« 

HH • -(E 4 AnR3*Fn»tSTn 
07 m (E -2.fAOR3*Fn*CTn 

BO ■ SORT( (0H-BT01**2 4 (0Z-HH|**2 4 BPn*»2» 

OH ■ BO - B 

IF( IPR.EO.l » WRITE(6.V5I(PDT( J.K»,J»l,5l,nH,r)Z.0iJ,Hn,H 
95 FORMAT!//. 2X.10F12.3I 


OuuuuaS' 
OLiOui'P'" 
0U*.*0 '*‘7 

OOOII'.*. W. 

IJIIU'IHC- ' 
Ul'l .< i« J ' 

(>1/ /DUHCI' 

0U0U'J95' 

HU'JOUOA 

'JUU’.'-P’ 

L VUC "97 ■ 

OiXXjiiPW" 

•jtni •!> ' > 
■10 vLi li'P 
(lO'l ') 

"JO ,Ml' ~ 
■K"H ' 7 

JlOX- 1 ' 

■«ij"j !o 7^ 

■ V I > 1 > 

I »t »r t, . 3 

‘M. « -1 I 
■ju ■ 1 1 i 
ui' ji n 

Ol,M.'Ol 1 .iv 
1 1 "Xl 

OOUOl I 5'.' 
0000 1 1 1^0 
000U117C 

uooullx.' 
0000 IIP. 
OW'iUl-!'" ■ 
Ji 

oo")ui2«;u 
oaoui/."* I 

M0')0i/9.| 
L'1 /uu 1 j; 5. 
Uil.jol/^ 
U00./I2 7. 
■lUOuUA J 
UU"Xj 129"' 




139 



sr.l ■ SCI ♦ 

i' .) 1 1 


Ml Til ino 


, II , ( •. . 

I" 

f 1 B hTCbsT ) 


U • ■ '1 


fy m MKBf,T5 

. 

UU. . i ' 


AfiTix m 


I'lt'i'i 1 A*. 


«L»>HA a Ft > 

* 

UJU il ?•» 


FI a At.FHA/^ 


"v‘JO 134 


F? a |K6 TF*AOHS»/k 


OUOC137 


sr.l a SC1 *P'*T<S.KI>iiF1 


tiu00l3fl 


sc? a S(.2 ♦ *JOT(S.K)aF7 


'J|) ||.il 3<- 


Sll » ill t»2*PI 


liUHK lAili 


S! <■ a SI? 4 F?BF? 


\ll;U 1 1 4l 


S*-l ■ SPl 4 Fl*Fl 


twJuiH,;i 

K'O 

CihmT I Nil® 


tiO 1 1 


Ji-( IS-i, g,i ) <;,i TD ?S 


1 / J 'IH., 


b a (Sr.jaSF2-Sr,?«‘Sll»/ISi!l*SI?-Sn®'«!n ) 



F| »-(*a<;Fi-«cn/SIl 


•■1 J(l 1 


xath' a Fi/fc 


UU.I Jl 47 


TTBIITH4 4( <CV-.;.IOS |a?A,n 


« 1 1 1 i '. ,1 • 4. V 


If,a|r.4l 


v»uU"l 


II^I IC.i.f.lSO » i;n TO »S 


U » h I * j * » 


«i<ITP(*',HOI 


.lUi'iiSl. 


FlK aT( lx. «=»<riM-. HfHg hAVF 

MOKE than 15C pr’I''iTS',//> 

7 >I 57'. 


!.*•( IC.c i.151 ) ICalSn 


• ■ •! A3 


Ii-( ic.t'i.ism c,n Tn son 


1,1, 1 ; 1 «H 

M’l 

bIMATI ir.llaTT 


iJl ■ i 


Fl'XATI Ii;,?)ae 


J. j 54 


ei«ATnr,.3»aFi 


•Ji.i .15 7., 


wKiTF<is,»»ni 


•'ij.il 5H.I 


FUKfAT(2X.»THIS IS Kjn.llTHH. TOTAL T I MF . OHM, LOCflL T.E.l*) 

uuCi.jl5Ks 


WkITF<A, 97» KOV.UTHH.TT.OftM.TML.f.F! 

OOOulft'JO 

<)7 

FnRMATI2X,I5.6PH.?» 


OOOUlFli.l 


IS-I a ? 


UU JU14?|. 


Stl a 0. 

i 

Ol'O'il a3i 


(ifi Tf' n 


ViO')014AJ 

2» 

A*- a U> 

' 

OUOOlASi) 


SCI ■ SHUT I sr.l /API 


tHJ'lul 441.1 


PkITEIA.^HI 


yO.JCl47.' 


F(mMAT(2X,*r)B>I.E, FI, RATIO. 

km SIGHAM 

!)Ui)0l4B0 


WRITE (i^.VAl OHM, E, FI. RATIO, SCI 

tiuajl4.li- 


FORMAT (5X,5FI2.5) 


Oil' 1.117 ■ > 


ISW a 1 


ULiUvl71' 


IPASS a INIH.n 


■K)u'il7<;i 


SCI a n. 


yuDiil 7.1 


SC2 » 0, 


-J0')ni74i. 
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sn » n. 

SI? ■ '1. 

S'*! ■ 

m 10(1. 

npM • ?M, 

|P « 0 

10 nwtAT ■ iNiv. n 
opi.iT»rwi,AT/I(i'). 

4r*> » AKSIOPLATI 

riH ■ Oh/10. 

I#lflOp.oT.2o.» (;n TO 40 
KsIaop.RT.OPf » r,0 TO ?0 
Ohm ■ 4i'P 

A( Sav ■ PM 5 .) i/ion 
ip(oh.(;t.ohm I r,n TO 4o 

11 OhI- ■ Orf 

AI.T • I •(4, I >41. 

OhL ■ OPtAT 
TM. . |m( 7. 1 I 

T( L ■ T'a/lOOO. 

«S ■ '(.SI >'11 » 

<nv«v ,(0 

4(1 lu ■ |tJ ♦ I 

IPnH.OT.lSni CALL AolMOlAol 
HT » 1(4,11 
ALdr. ■ lr;l3,l 1 
ALAT ■ IM(2, I I 
POT (I.IPI > 90. - DPLAT 
POT(?,IPI m HT/lO. 
poTn.iHi . OH 
P0T(4.IP1 > A'.AT/IOO. 

POT(S.IP| ■ ALOW/IOO. 

IP( IT(;ST.MS.O) r,n TO 410 

whits ( 4, 400) ALT.P0T(2, IPI.POT(3,1P|,POT(4, iPl.POTIS, IP) 
-400 POHMATIIX. 'ALT* SElO.a*' HT»*, 610.3, • 0«» ',610.3,' LAT»'. 
1610.3,' LONR-', 610.3) 

IT6ST-IT6ST41 
410 COnTIMOE 
I MTPaO 

XtAT»P0TI4, IP) 

HMf,HT«P0T<2,IP) 

XlOiitr.«AMnO(PDT(5,IP)43(SO.,360. ) 
IP(XLAT.L6.-50..0R.XLAT.G6.S0. ) 60 TO 440 
CALL ECOMPIXLAT.XLONG. HEIGHT, IHTP.OMAl 


'.•.I "") 7t 
1 1 7. 

U<» ■ *1 77 
r.i)..i(>l 7 A 
7“ 

• (OOOIA I 
00<JOl»l> 

B/ 

(HI H'lu J' 

• 44 
UlKIVl Ab' 
(.»' Jul 44 
.'C i-.l A7 

..,.J ««' 

u»* >' > I <(1' '• 

'■<’ (' 14| 

111 ) J 1 Hz’ 

■ i >1 0.7 
I 111 44 
I' ■ r ) U7 

f )-J*l 1 

■J.. ■■ 1I..7 
l•.••l',llu I 

OtiKjit 1 

il'J I'U i I. 

/ !' 

UU'iiJ/ I/' 

00002031 

(U).)'I204. 

PUi^(>2"3' 

'l(l■lu,:!l4l 

0UUP4>I7> 

0(|0J/'|«< 

UUUU,(UO' 

OOJ 0211 .H 

OOOU21K 

1HIJU21/.. 
OtH/U/l 3>. 
OIJi)U214. 

00002 161 ' 

UOJu^l 4t. 

000u2l7( 

Uv)Oii< 17^ 
UUUU21M0 
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•JO J'J ( t *■ 

c 


nw* |!i A.'MiMotv f.o««pr.Tlf^i 

• 

(tout / i 

c 


'.it « nM - OHA IN PfiTM,l*J| 


ti'l'iJi t ! 

f. 




Ui)')ti.f l'.' 



I<f|nn«.r,T,-'iO. |i»r>m.!t'i»>'OTIl» H»i-n**4 

. 

I'W 



C'lNTlNtlfe 





IP! ITfST.MI.il r,it TO htn 

- 

UU loi/** 



NKiTPKS.fcnl iKflTd. |p| 


OV't' //■ s 


h> » 

ITfST.ITfST+l 


OUOWXB'St' 


J 

PM«^4TC<IC, 'PlKol Hh |<; »,nn,in 




CiiiwT 


uOuotii- 



An Tn ‘i 


UOu"<ix- 


SOO 

W-ITMiS.4,Sni 


(■OUO.«i . ' 


hHi) 

P(iK.niT</,5*. 'ilMATn. Jl AS AllfSSf SPIS ITS/) 


gyou^ 3) ■ 



NX I TB ( *,, 7fiA I H f 1 <i«T ( 1 . .1 s 1 . 3 1 . 1 ■ 1 . r r. 1 


} // .•// 


7im1 

hUrt.'ATI 1 X,3f U.’l I 


UOUO<!43 



CA'.I. f!!'*SS3( PI> AT. ir.,ni)T2. If .4V») 
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Appendix Ci Dace reduction programs for estimating Q(o)) 

Includes t 

Had 

INt3 

Stack 

Imsl and system library functions are also required. 

Unlisted subroutines in Mad are listed In thesis by (Thayer, 1975) 


rtitt iAD 


rOITIAI A1 


•** Brona OnlvcrBitf Coiputar c«d* 


C lASIBTIC AIKLTSZS nOORiM fl' 

C OCTOIBR 1979 H’ 

C " 

DZRRRSZOR PRRZOD (25) ,0(500) ,RZRDOi (500) ,TZTLR(15) ,C(2) 

DZRBM8Z0I CORSO(2S) ,C0R(25) ,QRA81 (25) ,QRA62 (25) Rf 

OZRBNSZOI QPRAS 1(25) ,091882(25) B' 

C0R91.BZ Cl (1000), C2 (1000) ,9(2,2,25) ,110(1000) ,120(1000) R 

C0B9LIZ Ol(25),02(2S),tlRR8T(500) B 

DATA RAZB/1000/,RAZ9BR/25/ R' 

RADDBa-57. 295700 R' 

1 BRAD (5,1000,RR0>111) ZB90T,R9RR,ZRAI,l.ADI>,BBZR,8RL,RnRRT R/ 

1000 90RRAT (SZ5, 2910.0) R'' 

BRAD (5,2000) (9RRZ0D(Z) ,Z»1,R9RR) B< 

2000 PORRAT (0910.0) 

WRZTR (6,3000) R' 

3000 90RRAT (1R1,13(**») /I*,**' »11*»**VU,** RA88AT »*/ X, B/ 

1 •••,111,'*V1X,13(***)////1Z, 

2 'THR ZR90T DATA 8IRZR8;V/) R' 

c R? 

C P9T C0R99ZCZRRT8 ARB BRAD ZB R' 

c r; 

CALL Z099T (ZH9aT,rZTLR,B,RC0R9P,RAZR,DT,C1,C(1) ) B! 

CALL Z099T (ZR90T,ZZTLR,R,NC0R99,HAZR,0T,C2,C(2) ) R/ 

C BJ 

Z9 (ZRAZ.RO.O) ZRAZ>B R' 

RRZTB (6,11000) (C(Z) ,Z«1,2) ,8RL,LADD,ZHAZ,RRZR,R9BR, R/ 

1(9RRZ0D(Z) ,Z«1,R9RR) R' 

4000 90RBAT (///1Z,*TBR C0H90RRNT8 08RD ZB TRZ8 ARALmS ARE: X1>', R; 

1A4,' Z2-*,A«/////11, *68088X88 89RCTRAL RXBD0B8 : • //1Z, Bt 

2*8RLZCTXfXTT X8*,98.«/1Z, Bf 

3*B0BBRR or ZBTRR90LATX0R8 BBTBBRB BACH 9AZR*/1Z, Rf 

*•09 ADJACRBT *, BF 

«*RARHOBZCS Z8 89RCZ9ZRD A8* ,X3/1Z, RF 

5*THB BORBRR 09 RARB0RZC8 08RD OB BACB SXOB 09 RACB ZBTERPOLATZOR* , BF 
6* ZS*,X«/1Z, RF 

7*rRBRST TOTAL BRRR6T DRR8XTT ULORS ALLORBD ZB RACB SPECTRAL*, RF 
0* RZBDOW ZS*,ZB/1i:, BF 

9*TBBRB ARB*,Z3,* 9BRZ0DS ARALTZRD: */(1Z, 10912.2) ) RF 

TO«B*DT ‘ BF 

RRZTB (7,7000) (TZTLB(Z) ,Z*2, 15) RF 

7000 rORRAT (1«AR) RF 

c R* 

C TRARSXRRT SPBCTRA ZB TBB rRBQORRCZ BARDS ARB CALCOLATRD RF 

C BF. 

DO 10 L-1,R9RS RF 

RRZTB (6,17000) RF 

17000 rORRAT (1H1,72» ••) /1Z,72(**')/1X,20Z,***** BARD ABALTSZS ••••' RF 
1////) RF 

RRZTB (6,13000) (TZTLR(Z) ,Z-1,15) ,9BRZ0D(L) ,SRL BF 

13000 RORRAT (IB ,15A«//1Z,5Z, • PRRZOD** ,912. 2 ,* SBC.*,5Z, RF 

1*SBLBCTZfZTT«*,97.4) BF 

CALL 99ARHT (R,DT,BHXH,LAD0,ZBTBR9, 9BRX0D(L) ,SRL,R, RF 

1 XLOR,ZBZ6B,RrrT,BBBLOR,BAB0TB,RZTRRT) RF 

CALL HCOV (C1,C2,RRRR6T,XL0R,ZBX6B,B) RF 

C RF 
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mit IID rOKlIf A1 ••• Bcovb OBivarsitf coap«t«t Ctn* 

C IQOmOR {7|, PBOB 229 l.l. TUriR il.D. THBSXS, 1975 H 

C R 

RO-Zm fS.0132S6S*SBZ*R*DT/RIRZOD(l)»1.) R 

HO-ZPZZ (nOlT(RO) / R) H 

CUL RRRD (R,OT,ZR*Z,H,Q,ZI.OR«ZRZaR,RPPT»ZRBRR,RRRLOB,RAROrB, H 

1 C1,C2,Z1Q*Z2Q) H- 

C Bl 

C SHOOTIBD RRIRaZRS UR CUCOLATRD H' 

C Hf 

00«TO/RRRXOO(X,) R/ 

cm RRHOOR (Q,H,00«8RL«RXHDOR,9RXSRY1 Rf 

CALI PRRaRT (X1Q,Z1Q,RXH009,RRX8RT,B,P(1,1,L) ,RHRR8T) R' 

CALL FRR6RT (Z1Q,Z2Q,RXRDOB,RRXOHT, R, R( 1,2,L) ,RRRR6T) Rf 

p(2,i,Li>coHja (P(i,2,L)) r; 

CALL PRH8RT (Z20,Z2Q,RXIOOR,RRX8RT,R,r(2,2,L) «RHRR6X) R; 

c r; 

C CALCOLATR R8TXRATRS POR Q ■ X/R R! 

C R? 

Q1(L)-P(1t2*L)/P(1,1,L) H; 

Q2(L1«P(2,2,L)/P(2«1,L) Rf 

CORSQ(L)-RRAL( (P (2, 1,L) *P (2, 1,Ln /CP(1 *1* *P (2,2,L) ) ) R; 

COR (L) -SORT (CORSQ(M) R? 

QRA81(L)-SQRT( RRAL ( Q1 (L) *C0RJ8(Q1 (L)) ) ) Rr 

QPEAS1(L)>ATAR2( AZRA6(Q1(LH « RRALfOKL)) I * RADDR6 R/ 

QRA82(L)>SQRT( RRAL ( Q2 (L) *C0RJ8 (Q2 (L) ) ) ) R« 

QPRAS2 (L) -ATAM2 C AZRA8(Q2(LM r RRAL(Q2(LM ) * RADDR6 R; 

RRXTR (6,14000) RO,H Ht 

14000 PORRAT (///U,*TRRRR ARR DR8RRRS OP PRRR00H',3Z, Rf 

1* (R • SP8.4,') •) Rf 

RRXTR (6,14500) COH(L) Rf 

14500 PORRAT (/IX, 'CORRRRRart *,P12.4) Rf 

RRXTR (6, *.7000) Rf 

RRXTR (6,2005) RO Rf 

2005 PORRAT (////IX, *SROOTRRD SPRCTRAL RHRR6T RATRXZ*,5X,' ('*X4, Rf 

1 * DRORB'iS OP PRRROOR) •///26Z,*RZT*,27X,*XRTV9X,*.',60(* •)/ Hk 

29X,* •) j RA 

RRXTR (6,2001) V f1* 1*D «P(1»2*M Rk 

RRXTR (6,2002) P(2, 1,L) ,P(2,2,L) Rk 

2001 PORRAT (5X,'RXT *,2(2X,2R14.5) /9X, • •) Rf 

2002 PORRAT (5X,'XRT *,2(2X,2R14.5) /9X, • •) Rf 

H02-2*R0 RA 

CALL CORPXO (R02,C9RSQ(L),P(1,1,L),P(1,2,L),P(2,1,L),P(2,2,L) , Rf 

1 QRAOIf.L) ,QPRAS1(L) ,ORLQ,OBLP) Rk 

RRXTR (6,2004) Q1 (L) ,QRA61 (L) ,DBLQ,QPRAS1 (L) ,ORLP Rk 

2004 PORRAT (///IX,* RRSPORRB PORCTXOR RSTXHATRS t '//IX, Rk 

1'QI - ('.•2R14.5,«)'/20X,*RA6RXTODR - • ,P12. 4, IX, •«/- *,P12.4/20Z, Rf 

2 * PRASR • •,P12.2,1X,*t/> ',P5.2,' DR6RBES') Rk 

CALL CORPXD (R02,C0RS0(L) ,P(1,1,L) ,P(1,2,L) ,P(2,1,L) ,P(2,2,L) , Rk 

2 QRA62(L) ,QPRAS2(L) ,0RIQ,0RLP) Rk 

RRXTR (6,2006) Q2(L) ,QHA62(L) ,DRLQ,QPRAS2(L) ,DRLP Rk 

2006 PORRAT (1R0,'Q2 - (*,2R14.5,«) V20X,'RA6HXT0DR - *,P12.4, Rf 

1 IX,**/- *,P12.4/20X, Rf 

1 • PRASR ■*,P13.2,1Z,**/- *,P5.2,* DRORRRS*) Rf 

RRXTB (7,7001) PERIOD (L) ,SRL,RO,CORSQ(L) Rf 

RRXTB (7,7002) P(1, 1,L) ,P(1,2,L) Rk 
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PUBt BAD 


POmiB A1 


*** Bcevn OBivarsltf Coiputar Cant 



IBXTB (7,7002) P(2, 1,A) ,P(2,2,L) 
rOBBBT (P10.1,P10.2,X10,P10.5) 

h; 

7001 

n; 

7002 

POBBAT (2(2B,2B1A.5)) 

H' 

10 

COMTXBOB 

r; 


00 TO 1 

H' 

111 

STOP 



BBD 

R> 
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IMAO INT3 AONTMAN Cl TFST A/IA/AO 2II?A 

C Tf> APFt.V tiNfAH IMTfiRPOLATIDN TO MISSINR DATA POIHTS 1 I' 

C 1A.19A0 1 

C It. 

OIMiNSinN Tn2H».EIl?A>.MI12AI«CODE(12A> . ' <•' 

DIMENSinN ETn?"».FnU?A» "• 1 > 

oimenskin Ti«m.T?n5» , im- 

AEAO J9,inO» (Tim.I>ltl5l I 'T" 

AF AO 19*1001 IT2(||.|b1,1S| I>T<. 

RAAO 19.200) M I-'Tj 

ino POHMAT I15A4I 

2 nn AORIAT 1191 

RFAO (9.250) (T( I ).E( I ).FI 1 1 ),r.nOE( I ).E3(I).E19( I ). I>)«M) I 

?90 FfMM«T(lX*AE10.9) 

no 10 I-2.N !' T ■ 

IP )cnn*( D.Eo.o.) on to lo 

IF (COOEd-D.EO.l.) fio TO 10 1*1'. 

no 20 ,)»l.N 

IF (cnoEi l4.)).E'>.o. ) r,n to to I'Hf 

20 CONTIWIE 1'*>" 

30 AI-(E(I4.I)-E(I-1))/(T(I*J)-T(I-1)) IT. 

*9b(F 9( I4..l)-F9( 1-1 ) )/(TI I4J)-T(I-1 ) ) I I 

A?«(Fin4J)-FI( I-l ) )/)T(l4j)-T) I-l ) ) 1' I ■ 

AA-(FI9I 14.))-Fmi-1))/)T(I4.))-Tn-1)) I' " ' 

)»1«F11-1)-A1*T)1-1) • I*li- 

A9«E3(I-1)-A9*T( I-l) INTO 

R2»F1(I-1)-A2'*T(I-1) INT'J! 

»A-FI3)I-l)-AA*T(I-l) fwTO' 

on 40 K-i.j r'Ti!.' 

E( I-l4<)»Al*T( I-14K)4HI 1*1'. 

F3n-l4K)»A3t>T( I-14K)4«3 I' l'J 

Fl( I-l4<)»A2*T( I-14K)4H? 1''Tt‘ 

F13( 1-14K)»A4*T( I-14K)4«4 I' ft-. 

40 CONTINUE IMtU-j 

10 continue I 'Tu' 

WRITE (4.100) (Tl( I )*I«1.19) I'^T.Hj 

WRITE (4.100) (T2(I).I>1*19) , I^Tui) 

WRITE (4.200) N InTmu 

WRITE (4.300) (T(I)*E(I)*FI(I)*COOE(I)*E3(II.FI3(I).I>l.N) IM0(. 

300 FORMAT (AE10.3) l’'T' 

STOP 1^'T ) 

END 
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STACK FOKTKAN Cl TfST 1/2A/AO 9n<* 

stack 

nCT ?A,1979 

TO SK'.l-CTIV€t.V STAK EAiEKaies FO« A RIVEN PEKIOn OK 
FKEOUENCV AND WElfiHT EACH ENTRY BY ITS NO. Of OERKEfS OE PeeFiinf.. 
Nx«Nn. OF PEKinns CONSinEHED 
NSET-NO. OF 0AT4 SETS CONSIOERED 
NfiKMin-Nn. OF DATA SETS STACKED AT PEKIOD PEP(I) 
ClimtCi 2 m.C 2 im.C 22 m ABE A SEPARATE ENERRIV FIJNCT IOmS 
AT PEK( I » 

H.NO.OF DERREES OF FREEOOM FOR SPECTRAL ESTIMATES AT SnMc pfu( 
I I 

MO. OF OERHEES OF FREEOOM FOR PE«m 
NOSRI II«PINAt. NO, OF OERHEES OF FREOPM FOR PERm 
•w IF HISH A SkECTRAL element NOT TO HE INCLOnFO IM STACKINR 
AT THAT P6R« I » , SET N«n 
DIMENSION PER r 25 I . norm ( ?5 » . NDER ( 25 I 
DIMENSION SUMI(25>.SHM2(25I. S(JM3( 25 I . SUMAI 25 > 

DIMENSION SELI25I 

COMPLEX SCI U 25 > . SC 1 2 ( 25 I . SC21 ( 25 > . SC?2( 25 » 

COMPLEX Cll(25».C12(25I.C2l (25I.r,22(25> 

DATA SCll/25*«0.n,O.0>/,SCl2/25*(O.0.0«O)/,SC21/25'f(O.o,l/.i;)/, 

1 SC22/25«(0.0,0.0)/ 

n«TA St»Ml/25KO,0/,StlM2/25*0.0/,SUP3/25*0.0/.SOMA/25'»().0/ 

DATA NDER/25*0/,N0RM/25<'0/ 

NX»9 

NSBT-B 

RADnER«lBO,/3.1A15 

on 20 IX-l.NSET , 

DO 10 I-l.NX 

READI5.100) PERm.SELm.M 
100 FOKMATIFIO.I.EIO, 2.110) 

REAO(5.200) Clim.C12( I ) 

200 FnRMAT)2(2X.2ElA.5M 

RFAD( 5 . 200 ) C2im.C22m 
IF(N.NE.0)N0RM( 1 l-NORHI I )+l 

scnm>sciim4N*ciim 

SC12( 1 1-SC12I I l4N*C12n I 
sc2im>sc2im4NAC2i(i i 
Sr,22M i-SC22( 1 >4N*C22( 1 1 
NneRm«NDERm4N 
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1') CtiNTl'ine 

2n r.n.JTrnF 

nn ■^0 I ml, MX 

NhF(.( I lmwrvg(,( n*? 

w«nc(*«.250) •>«'«( 1 1 tNDfiRni.wBMin 
>‘5.1 FiiKiaTnx,»peHini)m»,Fl(t.1,' neRHtES OF F«FFnnf'»»,15.' STOH-'S 
w«iTc(*.,'?nn» H I) ,j.r.i 2 n i.sr, 2 i ( n.sc? 2 ( n 

^n^' FiwiAT(5x, 'cnmt.aeis.** 'Ciz-'tZEis.^./fSx.^cpi-'.^eis,*. 

1 •r,??m«,?=lS.A| 

fjl?mf:(mT(«Ffl!,(SCl2l I » I > »*“»?» 

m ?miU 2 /SOPT ( KgAL( sen ( 1 » >**2 > 

071m';(.'«T(.!FAUSr,22( I > »*•?! 

n?lmil?l /SUKT ( tiSAL ( SC?1 ( I) »**744 1 MAG « SC21 ( I I )**7> 
ph4sf«atam( AiPrtGtscizn n/K64Lisa2( i » ) iPKAnnEr, 
rflHm«SA..( ! 1 I A 1 WAG ( SC 1 2 I I) 

CfiHmC'iH/iWEALI 'iClK I I »'mFAL(SC22n ) I » 
w«ITF(A.4noi ni2,n?i .PMASe.r.nH 

Ano Ft.KV,flT( lx, •!*>» ..i(Cl?tm',F12,A,* 0(r,71 »m»,Fl2.A,» PHASE"'. F17, 
rnw7m •,=12.A» 

call criNFinff ngr,( | ) ,C0h,scil.Sr,12.SC71 .SC7?,oi?,PHASF.r>gL'i. 
••'t-lTi:(x.,x.noi ngt.o,r»ELP 

son FIIH. AT(5X, »ngLTA Ot. •,F12.A,' DEI TA PHim •,F1?,4./) 

'in c(i.MTrii!s 
STiiV 

•A" n 



Appendix D; Programs for calculating Q(w) 

Program name: Multi 

Subroutines: 

Solve 

Slmul 

Runge 

Step 

Safe 

Requires IMSL library and systems library 


1 • UIIAT (OJtOKIltIB. SYM'^) BHE XOCAL IW LB N17 S 8 72 


2 

C 

lJi)ITl 

3 

C 

JUNE 30,1900 

4 

C 

TO COHPUT MESrONSE FOR HULTl LAYER HODEL 

S 

C 


6 


rKC'GRAIi OCALUISP) 

7 


CUL CHANGE! ■4R0N‘‘> 

a 


CAIJ. At78icN<3,i8.'‘Oirrpirro''> 

9 


HEAL I1N,I2R 

10 


1)1 EHSJION X(4.4).UEFrHn(50) 

II 


IMi'lENHlON SlGnA(50).C<5e>,H(50),ARF<50),CNF<80) 

12 


N«3 

13 

C 


14 

C 

HEAD IN CONDUCriVITY HODEL 

10 


VariE<89,l) 

16 

1 

FOa'iAT<lX,“NO. OF SIGHAS?") 

17 


RE\rt(59,2) NS 

IB 

2 

FOJUIAT! 13) 

19 


DO 10 I«I,NS 

20 


WHri-ECOO.O) / 

21 

3 

FOHHATdX.'SIGHA AND DEPTH OF CHANCE IN HETERS*) 

22 


R£AD(59.100) SICMA( 1 ) .DEFTHHC 1 ) 

23 

100 

F0nnAT<2Fie.3) 

24 

10 

COmNUE 

23 


WHITE(U.4)(SIGnA(l).DEPTHN<I),l>l,NS) 

26 

4 

F0RHAT(8X, -SIGMA ".Eie.G," DEPTHH -.ElO.a) 

27 

C 


2B 


DO 500 JFR£0>1,19 

29 


X(1 ,I)>.1 

30 


X(2,l)>.2 

31 


X<3,l)>.3 

32 


X(4,l)>.4 

33 


X(1 ,2)>.5 

34 


X(2’,2)«.6 

30 


X<3,2>».7 

36 


X(4,2>«.8 

37 


X<l ,3)>.93 

38 


X(2.3)«.84 

39 


X(3,3)>.63 

40 


X(4,3)*.41 

41 


X(1.4)«.I6 

42 


X<2,4)».37 

43 


X(3,4)>.69 

44 


X(4.4)*.99 

43 


P« (FLOAT! JFRE0))»10. *44 

46 


ir(JFnEO.GE.ll) P>(FL0AT(JFRB0-9))»10.»«3 

47 


VRITE(.3.301 ) JFREO.P 

48 

501 

FORMAT! IX,'* JERBO M3,' ** INITIAL PERIOD -.EIO.S) 

49 


. EIN«I . 

30 


E2N«1. 

31 


GALL STEP!P,N,R,NU,C, IMAX, SIGMA, DEPTHH, NS, ZIN) 

32 


VRITE!3,III) IHAX,NH 

S3 

111 

FOnilAT!lX,'IMAX> ',13,' AND N0. OF STEPS « ',13,/) 

54 


CALL SAFE!C,NH,H.NH.N,ANF,CNF,IHAX,X,ZIN) 

33 

6 

FORMAT! IX, "USE ANF» *,E10.3," AND CNF* ',E10.3) 

56 


IlN»EIN/2.-ANF!IHAX) 

57 


I2NsE2N/2.-CNF(IMAX) 

58 

7 

F0RIIAT!1X,"I1N •,E10.3,"I2N ',E10.3) 

59 


Qsl 1N**2-H2N**2 

60 


0*SQRT!0/!E1W*E1H+E2N 'E2H) ) 
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61 l'flA8E«(IIN £2>^>I2N«<EIN) 

62 EHASE>ATAN(PHASr-^(I lN»EiN-«^12i<»E2N) ) 

63 WRITEO.n) (l,rKA!?E 

64 a FOllMATdX,"** /Q/ «",EI0.3,» PHASE* •.E10.3,/'/) 

68 500 COirriHUE 

06 CALL EXJT(I) 

07 END 
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6a 


SUBROUTINE SArElG, IG.H, IH.N.ANF.GNF . IHAX.X.ZIN) 

69 

C 

JULY 7,1980 

re 

G 

UFWARO INTEGRATION, 4 SOLUTIONS, SOLVE 

71 

C 

MATRIX FOR 4 GONSTANTS 

72 

G 


70 


DIIIENSION G( IG) ,HI IHI . ANF< IHAX) ,GNF( INAX) .Y(4) 

74 


DItlENSION X(4,4> 

70 


DIMENSION AN(80,4) ,CN(80,4) ,DAN(80,4) ,DGN(80,4) 

76 


DIMENSION DANn80), DGNF(80> ,A<4,4) ,8(4) ,VKAREA( 161 

77 

G 


78 

G 

** INITIALIZE 

79 


EIN*I. 

80 


E2N*I. 

8f 


IDCT»3 

82 


ALPHA*0.eOI 

83 


8ETA*0.00I 

84 


M>1 

88 


IA«4 

86 


IN«4 

87 


J«l / 

88 


V1UTE(3,100>N 

89 

100 

F0RMAT(1X,*F0R N« M8,* SKINDEPTHS* ,/) 

90 

8 

GONTINUE 

91 

no 

F0RMAT(1X,''READ IN INITIAL AN,DAN,GN,DGN *,/> 

92 


AN(I ,J)>X(1 ,J> 

93 


DANd ,J)«X(2,J) 

94 


Gild ,J)*X(3,J) 

98 


DCNd,J)>X(4,J) 

96 

120 

F0RMAT<F10.3) 

97 


Yd)»AN(I ,J) 

98 


Y(2)«DANd ,J) 

99 


Y(3)«GN(1 ,J) 

100 


Y(4)»DCNd ,J) 

101 


Z>ZIN 

192 

128 

F01UATdX,*Y(K) *,4E10.3) 

103 

G 


104 

G 

SOLVE SIMULTANEOUS EQUATIONS 

108 

G 


196 


CALL SinUL(Z,H,IH,N,G,lG,Y,ANF,GNF,DANF,DGNF,IHAX) 

107 


DO 10 1«1 .IMAX 

108 


ANd ,J)»ANF(I ) 

109 


1)AN(1 ,J)«DANF(1 ) 

no 


CNd ,J)«CNF(1) 

111 


DCNd ,J)*DCNF(1) 

112 

10 

CONTINUE 

113 

130 

FORMATdX,"! ",18," J *,I8," AN *,£10. 3,- DAN *,E10 

114 


1" CN ",E10.3," DCN *,£10.3) 

118 


J«J+I 

116 


IF<J.LT.8) GO TO 8 

117 

G 


118 

G 

DEFINE ELEMENTS OF MATRIX 

119 

G 


120 


DO 30 Jsl .4 

121 


Ad ,J)>ANd ,J) 

122 


A(2,J)«CNd ,J) 

123 


A<3.J)BDANdMAX,J)->^CdG>«ANdNAX.J) 

124 


A(4,J)«DGN(IMAX,J>-^GdG)»CNdMAX,J) 

123 

30 

CONTI NUE 

126 

200 

FOmiATdX,/,- ELEMENTS OF MATRIX*,/) 

127 

208 

FORMATdX, 4E10. 3) 


/,5X, 
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I2B 


129 


130 


131 


132 

210 

133 


134 

220 

135 


136 

230 

137 

C 

138 

C 

139 

C 

140 

C 

141 


142 


143 


144 


145 


146 

40 

147 

C 

140 

C 

149 

C 

130 


151 

800 

152 


153 

310 

154 


155 


156 


157 

320 

ISO 


159 

330 

160 


161 



D(l )«AU>I1A 
U<2)>DETA 

B(a)-0./2.»G<IC>«ElN 
B<4)>0./'2.i«C<lC)»E2N 
FOilHATUX," B<n> *.4EIB.9,/> 

CALL LBOTIF^A.n.lN.U.B.IDCT.VICAIlEA.IER) 

FOiU1AT(IX,*THE FOUR C0RSTANT8 ARE# *.4EiB.9> 

VlilTE(3,2aO> lER 

FOimATdX,/,*' ERROR PARAMETER ",I5) 

USING 1‘HEOREH THAT tIiERE EXISTS ONE SOLUTION , 

Yi «CI«YI l#C2»YI2#Ca)«iYI3#C4»YI4.ETC 

DO 40 1*1 .IHAX 

ANF(I )>B(1 >#AN(1 ,1)#B(2>»AN(1 ,2)#B(3)»ANn ,a)fB(4»XAN(l .4) 
DANFd )«B(i )»DAN(1 ,1>#B(2>»DAN(I .2)4^B<3>»0AN(l.a)#B(4)«DAN<l ,4) 
CNF( I )*Bd >»GN( I . I >4D(2)»CN(1 ,2)#B(3)«CN( 1 .3)4B(4)»GN(1 .4) 
DGNFd )«D< n«DGN(I . I )#B<2)»DCNd .2)#B(0)«DCNd ,3)#B(4)«DGN(1 .4) 
CONTINUE 

ARE COirriHUITY CONDITIONS MCT? 

WHITE<3,300) ALPHA, ANFd) 

FOnWATI/, IX. "ALPHA ".E10.3," ANFd) ",EI0.3) 
vniTEta.aio) beta.cnp(I) 

FORMAT (IX, "BETA ".EIO.3," CNFd ) ",E10.3) 

TESTI «(-l )#C( IC)»ANF< IMAX)+ (3./2.»C< IC)«EIN) 
TEST2*<-l)«CdC)«CNFdHAX)4'(3./2.»CdC)»E2N) 

WniTE(3,320) TESTl.DANFdHAX) 

FORMAT dX, "TEST I ", EIO.3. * COMPARED TO DAN< IHAX) ", EIO.3) 
V1HT£<3,330) TEST2,DCNF(IMAX) 

FOUMATdX, "TEST2 ", EIO.3," COfIPARED TO DCNdHAX) ", EIO.3) 

RETURN 

END 
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162 

16a 

104 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
ini 

IU2 
133 
184 
1 35 
186 

187 

188 
109 

190 

191 

192 

193 

194 

195 

196 

197 
190 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
21 1 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 


8UBItOUTI8E 8IHUL<Z.8,lH.N,CG,IG.y.Afl.GII.DA8.IICN.IIlAX) 

G ** AS OP APRIL 88,1980 

G CALCULATES SOLUTIOS TO PAIR OF COUPLED SIHULTAREOUS 

C 2ND ORDER DlPrCHEirriAL EOUATIONS OF FOUR D2A/DZ3* (8«C»»8 

C «AR/<Cw(Z-l>-fll»«84-2»ll«»a»CN. WHERE C IB A CIVEN COEF- 

C FlClEirr. RETURNS AN AND CN AS FUNCTION OP DEPTH. Z IS 
C NONDiriENSIONAL AND VARIES 0 TO I . REOUIRES FUNCTION RUNGE 
C 

DIRENSION AN<inAX>,CN(INAX>.Y(4),U(IH),F(4).CC<IC>,DAN(INAX> 
D I MENS I ON DCNdHAX) 

INTEGER RUNGE 
C 

C INITIALIZE 

C 

50 r0imAT(3X, **»» INITIALIZATION **•) 

60 FORnATUX.'Yl-.ElO.O,* Y2*.E10.3,* YS'.EIO.O.” Y4*,EI0.3) 
AN(I)«Y(I> 

CN<n«Y(3) 

DAN<1)>Y(8> . 

DCN(I)>Y<4> 

65 FORHATCIX," ANU )■ * ,EI0.3, * CH(I>- ■,E10.3./, 

65X,*DAN(I)> %EI0.3,'> DCN(1)» *,E10.3,/) 

I«1 

C 

C START INTEGRATION 

C 

4 IFIIR.EQi.l) STEP>HU) 

IFdH.EO.l ) GO TO 5 
5TEP>H<I > 

68 FORnATUX.*l -.IN,* STEP %E10.3) 

5 K>RUNCE(4,Y.F,Z.STEP) 

101 FORIUTCOX.-K'MS," Z>%E10.3) 

IFIK.EO.O) GO TO 10 

C 

C DERIVATIVES ARE DEFINED 

C 

IF(lC.Ea.l)C>GG(l) 

1F(IG.NE.1>C>CC(I) 

104 F0RHAT(1X,*C« *,EI0.3> . 

F(1 )»Y<2) 

F(2)*2«G«*2/(<G»(Z-1.)-1>1.)»«2»KY(I)+N«»2«2«Y(3) 

F(3)»Y(4) 

F(4)»2»C**2/< (C»(Z-l . )+l . >««2)»Y(3)-N»«2*2»Y< 1 ) 

102 FORMATnX, *'F1««,E10.3, " F2»»,EI0.3," F3«‘,E10.3,» F4«*,EI0.3> 
GO TO 5 

WHEN K IS RETURNED AS 0 .INTEGRATION IS COMPLETED 
10 1*1+1 

200 FORMATUX.-I* •,15,'* Z* *,E10.3) 

AN(I )»Y(1) 

CN( I )»Y(3) 

DAN(1 )*Y(2) 

DCN( I )»Y<4) 

300 F0RMAT(5X,"1»*.I5.*' AN(I ) * .EI0.3, ■ GN( I) * ,EIO .3) 

1F(I .LT. IHAX) GO TO 4 
IF(I.NE.IHAX) WR1TE(3,400) 

400 F0RnAT(9X." INDEX IMAX IS NOl' CONSISTENT* > 

RETURN 

END 
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C 

C 

c 

c 

G 

C 

G 

G 

C 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 


1ft 


822 
828 
824 
828 
22 « 

827 
228 
884 
288 
281 
888 
888 
284 
288 
28f 
887 
238 
889 
248 
241 
848 
248 
244 
248 
244 

247 

248 

249 
280 

281 G 

282 G 

283 1 

284 

288 G 

286 G 

287 2 
280 

289 

260 22 
241 
262 
268 
244 G 
268 G 
266 8 

267 

268 33 

269 

270 

271 G 

272 C 

273 4 

274 

278 44 

276 

277 

278 

279 G 

280 G 
201 8 


rUNCriON RUMGElM.Y.r.X.I) 

HE ruiiGTiof luiiGB nrtovi he pouim-GiiBai mimie - 
KOTTA nmMNi win ionTA*8 ooEmoium to ihtmiati a 
fViTEM or 2 •iHULTAiiioui fiiifrr 080ER oRoiRAinr oirrotEirriAL 
BBUATiom ruMinrui/iiK, AGiioii omt otep or 

LBNGH ■ 18 18B 18DEPB8DEirr VARIABLTii X. 8U8JICT TO 
IHITIAL GOIfOlTIOUf YU>, f J>1 ,8, . . ,81 . BAOI PfJl, ODll- 
VATIVE or YlJl. HUiT 8B COMPUTED POUR TlHEf.PDt HTTE- 
GRATI08 8TBP BY THE 0ALL180 PROORAM. THE rUHCTIOH HU8T BE 
CALLED PIVE TlNEi PER 8TEP IPAR81 1 ) . . .PA88 (8> > 80 TEAT 1VE 
I8DEPE8DE8T VARIA8LE VALUE X A8D TEE ROLUTIOR VALUIE 
<Y<1)...Y<8>> GAM RE UPDATED U818G TEE RUMCE-EUTTA AL- 
GORITHM. N 18 TEE PA88 COUHTER.RUMGE RBniR88 A8 1T8 VALUE 
1 TO 81GHAL TEAT ALL DER1VATIVB8 (TEE P(J>> RE EVALUA- 
TED OR 8 TO 81G8AL THAT TEE IRTERGRATIOH PR0CB88 PGR THE 
CURREHT 8TEP 18 PI818RED. 8AVEY(J) 18 U8ED TO 8AVE TIE 
THE IRITIAL VALUE OP Y(J> AMD PIICJl 18 TEE IHGREHEirT 
rUHGTKHI POR TEE J(TE) I0UATI08. A8 VR1TTE8. 8 HAY RE 
80 LARGER TEAM 88. CUUIRARAM.lJinEBR, ♦ VI 1XE8, APPLIED 
8UHER1CAL METIOINi. JOHN VILEY 4 8088, 8. Y., 1969. 

INTEGER RUNGE 

D1HEM810N rRl(08),8AVEY(88>,Y(8>,P(R> 

DATA M/8/ 

H*H41 

P0RMAT(7X. *PA88 MSI 

GO TO (1.8,8,4,8),H 


PA8S1 


RUNGE- 1 
RETTURM 

PA882 

DO 22 J-l.M 
8AVEy(J)-Y<J> 

PBI<J>-rCJ) 

Y < J ) - 6AVEV ( J ) 48 . 8»H4P ( J > 

X-X48.09H 

RUNGE- 1 

RETURN 

PA8S3 

DO 38 J-i ,N 

PHI ( J>-PRI ( J>42.8«P(J> 

y < J > -6AVEY < J >48 . 8«H»P ( J > 

RUNGE- 1 

RETURN 

PA8S4 

DO 44 J-l.N 

PH1<J>-PH1<J>42.09P(J> 

Y(J>-8AVEy(J)4H«P(J) 

X-X48.84H 

RUNGE- 1 

RETURN 

PASS8,.... 

DO 88 J-l.N 




k<)2 S8 Y(J)«SAVCy(J)4^(PHnJ>-^rCJ))»U/6.* 

2U3 H>0 

234 RUNGB»0 

208 RETURN 

206 C 
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288 SUBROUTINE STCP<P.N,B.NH.G,NZ,81GHA.DEPTilH,NSS.ZIN> 

289 C 

290 C SEPT. 4. 1980 (NOV28,|980> 

291 C TO DEPiNE STEP SIZES GIVEN CONDUCTIVITY PROFILE 

292 G P«PERIOD<SEG).N>SKIN DEPTHS NORMALIZATION. HU >> 

293 C STEP SIZE.NH-NO. OF STEPS, CC( I )>N»DELTA/(«371» 

294 C iO«»3).NZ«NO.OP DEPTHS Z>0 TO I.I.E. INAX IN SUB. 

298 C SAFE. 

296 C 

297 DIMENSION H(80) ,HH(80) ,DEPTH(80> .DELTAISO) .GC(80) .Z(80) 

298 DIMENSION SIGHA< 80) ,G(80> ,DEPTHM(50) 

299 NS«NS6 

300 PI«3.I4IS 

301 FREa«l./P 

302 DO 20 I«I ,NS 

303 DELTA! 1 )>SQRT(2.«I0.»»7/(0.«PI«»2»FRE0»SIGMA(I))) 

304 G(|>>N»DELTA(I)/(637I.»10.»»3) / 

308 DEPTHU)«DF.PTHM(1) 

306 IFU.EO.I) DEPTH(I)>l.>DEPTB(l)/(N»DELTA(m 

307 IFU.NE.I) DEPTHU)«DEPtH(I-l)-((D£PTHMU)-DEPTHM(l-l))/ 

308 I(N*DELTAU))) 

309 IFtCU ) .OE. I . ) DEPTH(I)«I. 

310 IF (DEPTH! 1) .GT. 1.) VRITE!89, 10) . I .DEPTH! I ) 

311 10 FORMAT! "EimOR I , 18, "DEPTH! I ) » *,£10.2) 

312 IF!DEPTH!I).LT.0.)NS«1 

313 IF!DEPTH(I).LT.O.) GO TO 21 

314 20 CONTINUE 

315 GO TO 22 

316 21 DEPTU(NS)*e. 

317 22 NH>0 

318 WRITE(3,200)!J,DEPTH!J),DELTA!J),C!J),J«I ,NS) 

319 200 FORMAT! IX, "J "U3," DEPTH ",E10.3,» DELTA ".ElO.O.'C *,EI0.3) 

320 DO 30 I«I,NS 

321 DC 25 J«1 ,5 

322 IFU.EO.I) HH!NH-«-J)*(l .~DEPTHU))/S. 

323 IFU.NE.I) UH!NH-^J)«!DEI'TH!l-l)-DEPTH!I))/8. 

324 CC(NH+J)»C!I ) 

325 IF(HH(NH+J) .EO.O) CC!NH-»^J)«C!I-U ) 

326 lF!HH(NH-t-J) .Ea.0)VRITE!3,5)!NH-)-.T),l ,DEPTHU) 

327 5 FORMAT! IX, "SPECIAL CONDITIONS NH4J«", 13," I •,I3,"DEPTH ".ElO.3) 

328 C »*NOTE CONDITIONS ON STEP HH ** 

329 25 CONTINUE 

330 NII«NH-«^5 

331 30 CONTINUE 

332 Z(I )«DEPTH(NS) 

333 300 FOlUIAT!/. IX, "THERE ARE ".lO," STEPS",/) 

334 DO 40 I«l ,NH 

335 K»NH+1-I 

336 M(I)«ini!K) 

337 C(1)*CC!K) 

3J8 ZU + l )=ZU)+H!I) 

339 400 FOUTIATUX,"! ".IS,*2« ",EI0.3."T0 Z* ".EIO.O," H= ",E10.3) 

310 40 CONTINUE 

341 NZ=mi+l 

342 IF<Z(NZ) .NE.I . ) !I(NH)»1 .-Z!NH) 

043 \>1UTE(3,450)(I ,ZU),Z!I-f| > ,H( I ) ,C!I ) , I » 1 ,NH) 

314 4 FOJUIAT!/, IX, "FINAL Z S ,H S, AND C S "./) 

315 450 F01UWT!IX."I ",I3,"Z« ",El0.3,"T0 Z= ",E10.3," H* ".E10.3, 

346 I " C' " .E10.3) 

317 ZIK=Z(I) 

318 RETURN 

319 f;nd- 
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390 

351 

392 

393 
354 
399 

356 

357 

358 

359 

360 

361 

362 

363 

364 
369 

366 

367 

368 

369 

370 
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G 

G 

G 

G 

G 

G 

G 

G 
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FUNCTION 

USAGE 

PARAMETERS 


C 

G 

G 

G 

C 

G 

C 

C 

C 

G 

C 

C 

C 

C 

C 

G 

C 

G 

C 

C 

G 

C 

C 

C 

C 

C 

G 

G 

C 

C 

C 

G 

C 

C 

C 


H 

N 

lA 

B 


IDOT 


WKAREA 

lER 


PRECISION 

HEQD. IHSL ROUTINES 
LANGUAGE 


LIBRARY 3 

- LINEAR B8UAT10N SOLUTION - PULL STORAGE 

NODE - SPACE ECONOMIZER SOLUTION. 

> CALL LEQTIF (A.N.N. IA.8. IDCT,VI:AREA, lER) 

- INPUT MATRIX OF DIMENSION N BY N CONTAINING 

THE COEFFCCIENT MATRIX OF THE EOUATION 
AX ■ B. 

ON OUTPUT, A IS REPLACED BY THE LU 

DECOMPOSITION OF A ROWWISE PERMUTATION OF 
A. 

- NUMBER OF RIGHT-HAND SIDES. ( INPUT! 

- ORDER or A AND NUMBER OF ROWS IN B. UNFIT) 

- NUMBER OF ROWS IN THE DIMENSION STATrurNT 

FOR A AND B IN THE CALLING PROGRAM. (INPUT) 

- INPUT MATRIX OF DIMENSION N BY M CONTAINING 

RIGHT-HAND SIDES OF THE EOUATION AX « B. 

ON OUTPUT. THE N BY H SOLUTION X REPLACES B. 

- INPUT OPTION. 

IF IDCT IS GREATER THAN 0, THE ELEMENTS OF 
A AND B ARE ASSUMED TO BE CORRECT TO lOGT 
DECIMAL DIGITS AND THE ROUTINE PERFORMS 
AN ACCURACY TEST. 

IF IDGT EOUALS ZERO, THE ACCURACY TEST IS 
BYPASSED. 

- WORK AREA OF DIMENSION GREATER THAN OR EQUAL 

TO N. 

- ERROR PARAMETER 
TERMINAL ERROR * 12B+N. 

N « 1 INDICATES THAT A IS ALGORITHMICALLY 
SINGULAR. (SEE THE CHAPTER L PRELUDE). 
WARNING ERROR « 32-)-N. 

N « 2 INDICATES THAT THE ACCURACY TEST 
FAILED. 

THE COMPUTED SOLUTION HAY BE IN ERROR 
BY MORE THAN CAN BE ACCOUNTED FOR BY 
THE UNCERTAINTY OF THE DATA. 

THIS WARNING CAN BE PRODUCED ONLY IF 
IDCT IS GREATER THAN 0 ON INPUT. 

SEE CHAPTER L PRELUDE FOR FURTHER 
DISCUSSION. 

- SINGLE 

- LUDATF.LUELMF.UERTST 

- FORTRAN 
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402 CC+ 

403 CC+ 

404 CC+ 

405 CC+ 

406 CC+ 

407 CC+ 
403 CC+ 
409 CC+ 
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DATE LAST 


ROUTINE: 

EDITION: 

CHANCED: 


LEQTIF 
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76 - 04-01 


CLASS TWO 


ROUTINES ARE MADE AVAILABLE BY NMC AS A SERVICE TO TI{E 
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a 


B 


C 

n 



c 


"i 

D 

DS 


Dst 

E 

E 




2n 


m 


en(t) 


f 


F( ) 
H 


^In’ ^2n 
l“(t) 


Glossary of Symbols 

radius of the earth 
series coefficient for p 

magnetic Induction 
series coefficient for p 

cross power spectra for series 1 and J 

cross power estimate 

normalized skin depth ■> N {/a 
constant coefficients 

electric displacement field; magnetic disturbance field 

disturbance local-time Inequality 

storm time variation 

activation energy 

electric field 

series coefficients for 6j^*e^(t) 

external magnetic field variation of order m, degree I 

frequency, cycles/sec 
differential equations 

fourier transform of 

magnetic field intensity 

series coefficients for ij^=i^(t) 

internal field variation of order m degree I 

free current density 
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M. 

N 

n, 



t, t 


n 


*^2j 


W(w,-w ) 
1 n 


X, X 

Y. 


count in j, index 

counting index; number of skin depths 
number of degrees of freedom 

scalar part of pololdal field S 

coefficients 

legendre polynomials 

response function 

estimate of response function 

electric charge 

radius vector 

radius 

selectivity 

pololdal vector field 

quiet time daily variations 

total time; total thickness 

time 

thickness of 1th slab (jth slab) with conductivity 

potential 
Gaussian window 

variable 

derivatives for finite difference algorithm 


z 


normalized radial distance 


Greek 


a 

8 

a 

Av 

Y- 


12 


'V 

e 


e 

X 

y 

0 

'V 

0 


P 

0 

Z 

e 


constant 

constant 

skin depth 

change In volume 

coherency between series 1 and 2 

arbitrary coefficient 

permittivity (electric) 
longitude 

permeability (magnetic); ion mobility 

phase 

function 

1/ay 

conductivity 

summation 

latitude 


0) 


frequency, 


rad/sec 
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